Cours 4 : « Prédiction et apprentissage automatique »

Sommaire

1. Prédiction par régression linéaire
2. (optionnel) Inférence dans le cadre de la régression

from IPython.display import Image
Image(filename='images/prédiction_et_apprentissage_automatique.jpg', width=500)

import pandas as pd
import seaborn as sns
import matplotlib
path_data = 'assets/data/'
matplotlib.use('Agg')
%matplotlib inline
import matplotlib.pyplot as plots
plots.style.use('fivethirtyeight')
import numpy as np

1. Prédiction par régression linéaire

Prédiction

Un aspect important de la science des données consiste à découvrir ce que les données peuvent nous dire sur l'avenir. Qu'est-ce que les données sur le climat et la pollution nous apprennent sur les températures dans quelques décennies ? Sur la base du profil internet d'une personne, quels sont les sites web susceptibles de l'intéresser ? Comment les antécédents médicaux d'un patient peuvent-ils être utilisés pour juger de sa réaction à un traitement ?

Pour répondre à ces questions, les spécialistes des données ont mis au point des méthodes permettant de faire des prédictions. Dans ce chapitre, nous étudierons l'une des méthodes les plus couramment utilisées pour prédire la valeur d'une variable en fonction de la valeur d'une autre variable.

Notre premier exemple est un ensemble de données comportant une ligne pour chaque chapitre du roman "Little Women". L'objectif est de prédire le nombre de caractères dans un chapitre (c'est-à-dire les lettres, les espaces, les signes de ponctuation, etc.) sur la base du nombre de phrases dans ce chapitre (dans la suite on utilise le nombre de points finaux comme une approximation simple du nombre de phrases).

little_women = pd.read_csv(path_data + 'little_women.csv')
little_women = little_women[['Periods'] + [col for col in little_women.columns if col != 'Periods']]
little_women.head(3)

	Periods	Characters
0	189	21759
1	188	22148
2	231	20558

sns.scatterplot(data=little_women, x='Periods', y='Characters')

<Axes: xlabel='Periods', ylabel='Characters'>

Comment pouvons-nous utiliser ces données pour prédire le nombre de caractères en fonction du nombre de phrases ?

Prédiction par régression linéaire

On peut imaginer plusieurs approches. Une des plus importantes à la fois historiquement et dans la pratique moderne en sciences des données repose sur modélisation linéaire des données.

L'idée est de prédire la variable d'intérêt $y$--appelée variable dépendante, ici il s'agit du nombre de caractères--comme une fonction linéaire, ou plus généralement affine, de la variable $x$ servant de base à la prédiction--appelée variable indépendante, ici il s'agit du nombre de phrases.

Notre prédiction $\hat{y}$ pour $y$ étant donné $x$ s'écrit alors: $\hat{y} = ax + b$, où $a$ et $b$, les paramètres du modèle, sont des coefficients à déterminer.

On appelle cette approche régression linéaire.

Etant donné un choix de valeurs pour la pente $a$ de la droiteet son intercept $b$ (aussi appelé ordonnée à l'origine), nous pouvons visualiser la qualité de la prédiction obtenue en superposant la droite de prédiction au diagramme de dispersion de l'échantillon utilisé. C'est ce que réalise la fonction ci-dessous.

def lw_errors(slope, intercept):
    sample = [[131, 14431], [231, 20558], [392, 40935], [157, 23524]]
    sns.scatterplot(data=little_women, x='Periods', y='Characters')
    xlims = np.array([50, 450])
    sns.lineplot(x=xlims, y=slope * xlims + intercept, color='blue', lw=2)
    for x, y in sample:
        plots.plot([x, x], [y, slope * x + intercept], color='r', lw=2)

Par exemple pour une pente de 87 caractères par phrase et un intercept de 4745 caractères nous obtenons le résultat ci-dessous. (On a fait figurer en rouge l'écart entre la valeur prédite et la valeur observée réellement pour quatre points de l'échantillon de données choisis arbitrairement).

slope = 50
intercept = 3000
print('Slope of Regression Line:    ', slope, 'characters per period')
print('Intercept of Regression Line:', intercept, 'characters')
lw_errors(slope, intercept)

Slope of Regression Line:     50 characters per period
Intercept of Regression Line: 3000 characters

La ligne ci-dessus sous-estime largement le nombre de caractères en étant située trop "bas." Peut-on trouver un meilleur choix de paramètre pour modéliser notre échantillon de données ?

Après quelques tatônnements on arrive à un meilleur résultat. Par exemple :

slope = 83
intercept = 4500
print('Slope of Regression Line:    ', slope, 'characters per period')
print('Intercept of Regression Line:', intercept, 'characters')
lw_errors(slope, intercept)

Slope of Regression Line:     83 characters per period
Intercept of Regression Line: 4500 characters

On peut lire sur ce graphe la prédiction du nombre de caractères pour n'importe quel nombre de phrases, y compris un nombre de phrases pour lequel on n'avait pas de données dans l'échantillon.

On peut aussi calculer la prédiction directement avec la formule $\hat{y}=ax+b$.

Le choix ci-dessus semble donner un bien meilleur résultat (une meilleure adéquation du modèle aux données), cependant la procédure consistant à choisir une bonne valeur des paramètres "à l'oeil" comme nous venons de le faire, n'est pas très satisfaisante.

Nous allons voir une approche plus systématique et très utilisée en sciences des données : la méthodes des moindres carrés.

La méthode des moindres carrés

Pour trouver les meilleurs paramètres, nous avons besoin d'une définition raisonnable du terme "meilleur". Rappelons que le but de la droite est de prédire ou estimer les valeurs de $y$, étant donné les valeurs de $x$. Les estimations ne sont généralement pas parfaites. Chacune d'entre elles s'écarte de la valeur réelle d'une erreur. Un critère raisonnable pour qu'une droite soit la "meilleure" est qu'elle ait l'erreur globale la plus petite possible parmi toutes les droites.

Dans cette section, nous allons préciser ce critère et voir si nous pouvons identifier la meilleure droite selon ce critère.

À chaque point du diagramme de dispersion correspond une erreur de prédiction calculée comme la valeur réelle moins la valeur prédite. Il s'agit de la distance verticale entre le point et la ligne, avec un signe négatif si le point est en dessous de la ligne. Sur le graphe ci-dessous cette erreur est représenté en rouge pour quatre points choisis arbitrairement.

slope = 83
intercept = 4500
print('Slope of Regression Line:    ', slope, 'characters per period')
print('Intercept of Regression Line:', intercept, 'characters')
lw_errors(slope, intercept)

Slope of Regression Line:     83 characters per period
Intercept of Regression Line: 4500 characters

Si nous avions utilisé une autre ligne pour créer nos estimations, les erreurs auraient été différentes. Le graphique ci-dessous montre l'ampleur des erreurs si nous avions utilisé une autre droite pour l'estimation. Le deuxième graphique montre les erreurs importantes obtenues en utilisant une ligne qui est carrément ridicule.

lw_errors(50, 10000)

lw_errors(-100, 50000)

Erreur quadratique moyenne

Nous avons maintenant besoin d'une mesure globale de la taille approximative des erreurs. Vous reconnaîtrez l'approche utilisée pour créer cette mesure - c'est similaire à la façon dont on peut expliquer la formule pour l'écart-type.

Si vous utilisez une ligne arbitraire pour calculer vos estimations, certaines de vos erreurs seront probablement positives et d'autres négatives. Pour éviter toute annulation lors de la mesure de la taille approximative des erreurs, nous prendrons la moyenne des erreurs au carré plutôt que la moyenne des erreurs elles-mêmes.

L'erreur quadratique moyenne d'estimation est une mesure de l'importance approximative des erreurs quadratiques, mais comme nous l'avons noté précédemment, ses unités sont difficiles à interpréter. En prenant la racine carrée, on obtient l'erreur quadratique moyenne (rmse), qui est exprimée dans les mêmes unités que la variable prédite et qui est donc beaucoup plus facile à comprendre.

Minimiser l'erreur quadratique moyenne

Les observations que nous avons faites jusqu'à présent peuvent être résumées comme suit.

• Pour obtenir des estimations de $y$ basées sur $x$, vous pouvez utiliser n'importe quelle ligne.
• Chaque ligne a une erreur quadratique moyenne d'estimation.
• Les "meilleures" lignes ont des erreurs plus faibles.

Existe-t-il une "meilleure" droite ? En d'autres termes, existe-t-il une droite qui minimise l'erreur quadratique moyenne parmi toutes les droites ?

Pour répondre à cette question, nous commencerons par définir une fonction lw_rmse qui calcule l'erreur quadratique moyenne de n'importe quelle ligne du diagramme de dispersion. La fonction prend la pente et l'ordonnée à l'origine (dans cet ordre) comme arguments.

def lw_rmse(slope, intercept):
    lw_errors(slope, intercept)
    x = little_women['Periods']
    y = little_women['Characters']
    fitted = slope * x + intercept
    mse = np.mean((y - fitted) ** 2)
    return mse**0.5

rmse = lw_rmse(50, 10000)
print("Root mean squared error:", rmse)

Root mean squared error: 4322.167831766537

rmse = lw_rmse(-100, 50000)
print("Root mean squared error:", rmse)

Root mean squared error: 16710.11983735375

Les mauvaises lignes ont des valeurs élevées de rmse, comme prévu. Mais la rmse est beaucoup plus faible si nous choisissons une pente et une ordonnée à l'origine plus raisonnable.

rmse = lw_rmse(90, 4000)
print("Root mean squared error:", rmse)

Root mean squared error: 2715.5391063834586

Nous avons maintenant un critère clair pour choisir les paramètres. Nous allons maintenant utiliser python pour chercher systématiquement une bonne valeur pour les paramètres. Pour ce faire, nous allons utiliser la fonction scipy.minimize du module scipy de python.

Optimisation numérique

La fonction minimize peut être utilisée pour trouver les arguments d'une fonction pour lesquels la fonction renvoie sa valeur minimale. Il y aurait beaucoup à dire sur le fonctionnement de cette méthode, mais on se contentera ici de dire qu'elle part d'une valeur de départ pour les paramètres et cherche ensuite comment légèrement modifier ces paramètres pour les améliorer. Elle obtient ainsi une nouvelle valeur des paramètres qui donne de meilleurs résultats. Elle réitère ensuite cette approche, autant de fois qu'il le faut, jusqu'à atteindre un point ou il n'existe plus de façon de modifier légèrement les paramètres pour réduire l'erreur quadratique.

La figure ci-dessous permet de visualiser le genre de point auquel on peut arriver par cette méthode. Il faut imaginer que l'axe des $x$ représente le choix du paramètre (dans cette visualisation il n'y a qu'un seul paramètre à choisir) et l'axe des $y$ représente l'erreur quadratique moyenne.

La procédure de recherche décrite informellement ci-dessus devrait nous mener à un minimum pour l'erreur quadratique moyenne, ou du moins un minimum local. Il n'est pas garanti qu'elle trouve un minimum global cependant, car elle pourrait se retrouver coincée dans une "cuvette" et rester au niveau d'un minimum local. On qualifie donc la méthode de recherche utilisée de méthode de recherche locale.

Pour des problèmes simples comme celui considéré ici, on pourrait utiliser d'autres approches (dans une section plus bas on trouve une formule explicite qui donne directement les meilleurs paramètres par exemple), mais dans les applications modernes des sciences des données (par exemple les grands modèles de language, les générateurs d'images etc.) les seules solutions utilisables en pratique sont basée sur des méthodes de recherche locales comme celle qu'on a considéré dans cette section.

La distinction entre minimum local et global et le problème de trouver des procédures de recherche locale qui évitent les minimum locaux sont donc très important. Pour le cas de la régression linéaire au moindres carrés qu'on considère ici, cependant il est possible de montrer que la fonction d'erreur considérée (erreur quadratique moyenne) possède un unique minimum. Il n'y a donc pas de risque de se tromper : si on trouve un minimum local, c'est forcément l'unique minimum global de la fonction d'erreur!

L'argument de minimize est une fonction qui prend elle-même des arguments numériques et renvoie une valeur numérique. Par exemple, la fonction lw_rmse prend une pente et une interception numériques comme arguments et renvoie la rmse correspondante.

L'appel minimize(lw_rmse) renvoie un tableau composé de la pente et de l'ordonnée à l'origine qui minimisent la rmse. Ces valeurs minimisantes sont d'excellentes approximations obtenues par essais et erreurs intelligents, et non des valeurs exactes basées sur des formules.

from scipy.optimize import minimize

def lw_rmse_no_plot(slope, intercept):
    # same function as before except that we don't make a plot
    x = little_women['Periods']
    y = little_women['Characters']
    fitted = slope * x + intercept
    mse = np.mean((y - fitted) ** 2)
    return mse**0.5

def lw_rmse_for_minimize(params):
    slope, intercept = params
    return lw_rmse_no_plot(slope, intercept)

result = minimize(lw_rmse_for_minimize, x0=[0, 0])  # Initial guess for slope and intercept
best = result.x  # Extracting the best slope and intercept
best

array([  86.97640358, 4745.11706869])

Ces valeurs sont les mêmes que celles que nous avons calculées précédemment en utilisant les fonctions slope et intercept. Nous constatons de petites déviations dues à la nature inexacte de minimize, mais les valeurs sont essentiellement les mêmes.

rmse = lw_rmse(best[0], best[1])
print("Root mean squared error:", rmse)

Root mean squared error: 2701.6907879382966

On obtient une erreur quadratique moyenne d'environ $2702$, qui est bien la plus faible qu'on ait obtenu.

print("slope from minimize:       ", best[0])
print("intercept from minimize:   ", best[1])

slope from minimize:        86.97640358176976
intercept from minimize:    4745.117068692665

Corrélation

La méthode de recherche de paramètre que nous venons de voir est très générale et très utilisée en pratique. Dans le cas particulier de la régression linéaire, il se trouve qu'il est possible de trouver une formule explicite pour les paramètres (pente et intercept) minimisant l'erreur quadratique moyenne. Notez que dans la plupart des cas en sciences des données il n'est pas possible d'obtenir de formules explicite et que les méthodes de recherche automatisée telle que celle considérée ci-dessus sont souvent les seules solutions envisageables.

Nous allons maintenant dériver cette formule explicite par une approche intuitive (sans démonstration) et, pour ce faire, nous commencons dans cette section, par développer une mesure du degré de regroupement d'un diagramme de dispersion autour d'une ligne droite. Formellement, cela s'appelle mesurer l'association linéaire ou corrélation entre deux variables (les variables correspondants aux axes du diagramme).

Le tableau hybrid contient des données sur les voitures hybrides vendues aux Etats-Unis de 1997 à 2013. Les données ont été adaptées à partir des archives de données en ligne du [Prof. Larry Winner] (http://www.stat.ufl.edu/%7Ewinner/) de l'Université de Floride. Les colonnes :

• véhicule : modèle de la voiture
• year : année de fabrication
• msrp : prix de détail suggéré par le fabricant en dollars de 2013
• acceleration : taux d'accélération en km par heure par seconde
• mpg : efficacité energétique en miles par gallon
• class : classe du modèle.

hybrid = pd.read_csv(path_data + 'hybrid.csv')

hybrid

	vehicle	year	msrp	acceleration	mpg	class
0	Prius (1st Gen)	1997	24509.74	7.46	41.26	Compact
1	Tino	2000	35354.97	8.20	54.10	Compact
2	Prius (2nd Gen)	2000	26832.25	7.97	45.23	Compact
3	Insight	2000	18936.41	9.52	53.00	Two Seater
4	Civic (1st Gen)	2001	25833.38	7.04	47.04	Compact
...	...	...	...	...	...	...
148	S400	2013	92350.00	13.89	21.00	Large
149	Prius Plug-in	2013	32000.00	9.17	50.00	Midsize
150	C-Max Energi Plug-in	2013	32950.00	11.76	43.00	Midsize
151	Fusion Energi Plug-in	2013	38700.00	11.76	43.00	Midsize
152	Chevrolet Volt	2013	39145.00	11.11	37.00	Compact

153 rows × 6 columns

Le graphique ci-dessous est un diagramme de dispersion donnant prix de la voiture msrp (en dollars) en fonction de sa capacité d'acceleration (en km/h/s). msrp est donc représenté sur l'axe vertical et acceleration sur l'axe horizontal.

sns.scatterplot(data=hybrid, x='acceleration', y='msrp')

<Axes: xlabel='acceleration', ylabel='msrp'>

Remarquez l'association positive. La dispersion des points est en pente ascendante, ce qui indique que les voitures ayant une plus grande accélération ont tendance à coûter plus cher, en moyenne ; de manière équivalente, les voitures qui coûtent plus cher ont tendance à avoir une plus grande accélération en moyenne. (Que pensez-vous de la direction causale de cette association ? Est-ce qu'un prix élevé cause une capacité d'accéleration élevée ? Est-ce qu'une capacité d'accélération élevée cause un prix élevé ? Aucun des deux ?)

Le diagramme de dispersion du prix msrp en fonction de l'efficacité énergétique mpg (le nombre de miles qu'il est possible de parcourir avec un gallon de carburant) montre une association négative. Les voitures hybrides ayant une meilleure efficacité énergétique ont tendance à coûter moins cher, en moyenne. Cela n'est pas si surprenant si l'on considère que les voitures qui peuvent accélèrer rapidement ont aussi tendance à être moins économes en carburant et donc à avoir un efficacité énergétique plus faible. Comme le montre le diagramme de dispersion précédent, ce sont également les voitures qui ont tendance à coûter plus cher.

sns.scatterplot(data=hybrid, x='mpg', y='msrp')

<Axes: xlabel='mpg', ylabel='msrp'>

Outre l'association négative, le diagramme de dispersion du prix par rapport à l'efficacité montre une relation non linéaire entre les deux variables. Les points semblent être regroupés autour d'une courbe et non d'une ligne droite.

Si nous limitons les données à la seule catégorie des SUV, l'association entre le prix et l'efficacité est toujours négative, mais la relation semble plus linéaire. La relation entre le prix et l'accélération des SUV montre également une tendance linéaire, mais avec une pente positive.

suv = hybrid[hybrid['class'] == 'SUV']
sns.scatterplot(data=suv, x='mpg', y='msrp')

<Axes: xlabel='mpg', ylabel='msrp'>

sns.scatterplot(data=suv, x='acceleration', y='msrp')

<Axes: xlabel='acceleration', ylabel='msrp'>

Nous pouvons donc tirer des informations utiles de l'orientation générale et de la forme d'un diagramme de dispersion.

Regardons à présent si l'orientation et la forme du nuage de points dépend des unités dans lesquelles les variables ont été mesurées.

Mesure en unités standard

Les unités standard mesurent le nombre d'écarts-types au-dessus ou en-dessous de la moyenne.

Certaines valeurs des unités standard sont négatives, correspondant à des valeurs initiales inférieures à la moyenne. D'autres valeurs d'unités standard sont positives. Mais quelle que soit la distribution de nos données, les bornes de Tchebychev impliquent que les unités standard se situent généralement dans l'intervalle (-5, 5).

Pour convertir une valeur en unités standard, il faut d'abord déterminer l'écart par rapport à la moyenne, puis comparer cet écart à l'écart-type.

$$ z ~=~ \frac{\mbox{valeur }-\mbox{moyenne}}{\mbox{SD}} $$

Comme nous le verrons, les unités standard sont fréquemment utilisées dans l'analyse des données. Il est donc utile de définir une fonction qui convertit un tableau de nombres en unités standard.

Définissons la fonction standard_units pour convertir un tableau de nombres en "unités standard", c'est à dire pour chaque valeur enlever la moyenne et diviser par l'écart-type.

def standard_units(any_numbers):
    "Convert any array of numbers to standard units."
    return (any_numbers - np.mean(any_numbers)) / np.std(any_numbers)

Nous pouvons utiliser cette fonction pour redessiner les deux diagrammes de dispersion pour les SUV, avec toutes les variables mesurées en unités standard.

standardized_data = pd.DataFrame({
    'mpg (standard units)': standard_units(suv['mpg']),
    'msrp (standard units)': standard_units(suv['msrp'])
})

sns.scatterplot(data=standardized_data, x='mpg (standard units)', y='msrp (standard units)')
plots.xlim(-3, 3)
plots.ylim(-3, 3)

(-3.0, 3.0)

standardized_data = pd.DataFrame({
    'acceleration (standard units)': standard_units(suv['acceleration']),
    'msrp (standard units)': standard_units(suv['msrp'])
})

sns.scatterplot(data=standardized_data, x='acceleration (standard units)', y='msrp (standard units)')
plots.xlim(-3, 3)
plots.ylim(-3, 3)

(-3.0, 3.0)

Les associations que nous voyons dans ces figures sont les mêmes que celles que nous avons vues précédemment. De plus, comme les deux diagrammes de dispersion sont maintenant dessinés exactement à la même échelle, nous pouvons voir que la relation linéaire dans le second diagramme est un peu plus floue que dans le premier.

Nous allons maintenant définir une mesure qui utilise des unités standard pour quantifier les types d'association que nous avons observés.

Le coefficient de corrélation

Le coefficient de corrélation mesure la force de la relation linéaire entre deux variables. Graphiquement, il mesure le degré de regroupement du diagramme de dispersion autour d'une ligne droite.

Le terme coefficient de corrélation n'est pas facile à prononcer, c'est pourquoi il est généralement abrégé en corrélation et désigné par $r$.

Voici quelques faits mathématiques concernant $r$ que nous observerons par simulation.

• Le coefficient de corrélation $r$ est un nombre compris entre $-1$ et 1.
• $r$ mesure la mesure dans laquelle le nuage de points se concentre autour d'une ligne droite.
• $r = 1$ si le diagramme de dispersion est une ligne droite parfaite inclinée vers le haut, et $r = -1$ si le diagramme de dispersion est une ligne droite parfaite inclinée vers le bas.

La fonction r_scatter prend une valeur de $r$ comme argument et simule un nuage de points avec une corrélation très proche de $r$. En raison du caractère aléatoire de la simulation, la corrélation n'est pas censée être exactement égale à $r$.

Appelez r_scatter plusieurs fois, avec différentes valeurs de $r$ comme argument, et voyez comment le nuage de points change.

Lorsque $r=1$, le nuage de points est parfaitement linéaire et s'incline vers le haut. Lorsque $r=-1$, le nuage de points est parfaitement linéaire et s'incline vers le bas. Lorsque $r=0$, le diagramme de dispersion est un nuage informe autour de l'axe horizontal, et les variables sont dites non corrélées.

def r_scatter(r):
    """
    Génère un nuage de points avec une corrélation approximative de r.
    
    Paramètre :
    r (float) : Coefficient de corrélation souhaité, entre -1 et 1.
    """
    # Génération de 1000 points pour x selon une distribution normale standard
    x = np.random.normal(0, 1, 1000)
    # Génération de 1000 points pour z selon une distribution normale standard
    z = np.random.normal(0, 1, 1000)
    # Calcul de y pour obtenir la corrélation souhaitée avec x
    y = r * x + np.sqrt(1 - r**2) * z
    # Création du DataFrame pour faciliter le tracé avec seaborn
    data = pd.DataFrame({'x': x, 'y': y})
    # Tracé du nuage de points
    sns.scatterplot(data=data, x='x', y='y')
    # Définition des limites des axes pour une meilleure visualisation
    plots.xlim(-4, 4)
    plots.ylim(-4, 4)

r_scatter(0.9)

r_scatter(0.25)

r_scatter(0)

r_scatter(-0.55)

Calculer $r$

La formule de $r$ n'est pas apparente d'après nos observations jusqu'à présent. Elle repose sur une base mathématique qui n'entre pas dans le cadre de ce cours. Cependant, comme vous le verrez, le calcul est simple et nous aide à comprendre plusieurs propriétés de $r$.

Formule pour $r$ :

$r$ est la moyenne des produits des deux variables, lorsque les deux variables sont mesurées en unités standard.

Voici les étapes du calcul. Nous allons appliquer ces étapes à un simple tableau de valeurs de $x$ et $y$.

x = np.arange(1, 7, 1)
y = np.array([2, 3, 1, 5, 2, 7])
t = pd.DataFrame({
    'x': x,
    'y': y
})
t

	x	y
0	1	2
1	2	3
2	3	1
3	4	5
4	5	2
5	6	7

D'après le diagramme de dispersion, nous nous attendons à ce que $r$ soit positif mais pas égal à 1.

sns.scatterplot(data=t, x='x', y='y', s=30, color='red')

<Axes: xlabel='x', ylabel='y'>

Étape 1. Convertir chaque variable en unités standard.

t_su = t.copy()
t_su['x (standard units)'] = standard_units(t['x'])
t_su['y (standard units)'] = standard_units(t['y'])
t_su

	x	y	x (standard units)	y (standard units)
0	1	2	-1.46385	-0.648886
1	2	3	-0.87831	-0.162221
2	3	1	-0.29277	-1.135550
3	4	5	0.29277	0.811107
4	5	2	0.87831	-0.648886
5	6	7	1.46385	1.784436

**Étape 2 : Multiplier chaque paire d'unités standard.

t_product = t_su.copy()
t_product['product of standard units'] = t_su['x (standard units)'] * t_su['y (standard units)']
t_product

	x	y	x (standard units)	y (standard units)	product of standard units
0	1	2	-1.46385	-0.648886	0.949871
1	2	3	-0.87831	-0.162221	0.142481
2	3	1	-0.29277	-1.135550	0.332455
3	4	5	0.29277	0.811107	0.237468
4	5	2	0.87831	-0.648886	-0.569923
5	6	7	1.46385	1.784436	2.612146

**Étape 3 : $r$ est la moyenne des produits calculés à l'étape 2.

# r is the average of the products of standard units

r = np.mean(t_product['product of standard units'])
r

np.float64(0.6174163971897709)

Comme prévu, $r$ est positif mais n'est pas égal à 1.

Propriétés de $r$

Le calcul montre que :

• $r$ est un nombre pur. Il n'a pas d'unités. En effet, $r$ est basé sur des unités standard.
• $r$ n'est pas affecté par un changement d'unité sur l'un ou l'autre des axes. Ceci est également dû au fait que $r$ est basé sur des unités standard.
• Le changement d'axe n'affecte pas $r$. Algébriquement, c'est parce que le produit des unités standard ne dépend pas de la variable appelée $x$ et de celle appelée $y$. Géométriquement, le changement d'axe reflète le nuage de points autour de la ligne $y=x$, mais ne modifie pas le degré de regroupement ni le signe de l'association.

sns.scatterplot(data=t, x='y', y='x', s=30, color='red')

<Axes: xlabel='y', ylabel='x'>

La fonction de corrélation

Nous allons calculer des corrélations de façon répétée, il est donc utile de définir une fonction qui les calcule en effectuant toutes les étapes décrites ci-dessus. Définissons une fonction corrélation qui prend un tableau et les étiquettes de deux colonnes du tableau. La fonction renvoie $r$, la moyenne des produits des valeurs de ces colonnes en unités standard.

def correlation(t, x, y):
    return np.mean(standard_units(t[x]) * standard_units(t[y]))

Appelons la fonction sur les colonnes x et y de t. La fonction renvoie la même réponse à la corrélation entre $x$ et $y$ que celle obtenue par l'application directe de la formule pour $r$.

correlation(t, 'x', 'y')

np.float64(0.6174163971897709)

Comme nous l'avons remarqué, l'ordre dans lequel les variables sont spécifiées n'a pas d'importance.

correlation(t, 'y', 'x')

np.float64(0.6174163971897709)

En appelant correlation sur les colonnes du tableau suv, on obtient la corrélation entre le prix et le kilométrage ainsi que la corrélation entre le prix et l'accélération.

correlation(suv, 'mpg', 'msrp')

np.float64(-0.6667143635709919)

correlation(suv, 'acceleration', 'msrp')

np.float64(0.48699799279959155)

Ces valeurs confirment ce que nous avions observé :

• Il existe une association négative entre le prix et l'efficacité, tandis que l'association entre le prix et l'accélération est positive.
• La relation linéaire entre le prix et l'accélération est un peu plus faible (corrélation d'environ 0,5) que celle entre le prix et le kilométrage (corrélation d'environ -0,67).

La corrélation est un concept simple et puissant, mais il est parfois mal utilisé. Avant d'utiliser $r$, il est important de savoir ce que la corrélation mesure et ne mesure pas.

Association n'est pas causalité

La corrélation ne mesure que l'association. La corrélation n'implique pas la causalité. Bien que la corrélation entre le poids et les aptitudes en mathématiques des enfants d'un district scolaire puisse être positive, cela ne signifie pas que le fait de faire des mathématiques alourdit les enfants ou que la prise de poids améliore leurs aptitudes en mathématiques. L'âge est une variable confondante : les enfants plus âgés sont à la fois plus lourds et meilleurs en mathématiques que les enfants plus jeunes, en moyenne.

La corrélation mesure l'association linéaire

La corrélation ne mesure qu'un seul type d'association - linéaire. Les variables qui ont une forte association non linéaire peuvent avoir une corrélation très faible. Voici un exemple de variables qui ont une relation quadratique parfaite $y = x^2$ mais dont la corrélation est égale à 0.

new_x = np.arange(-4, 4.1, 0.5)
nonlinear = pd.DataFrame({
    'x': new_x,
    'y': new_x**2
})

sns.scatterplot(data=nonlinear, x='x', y='y', s=30, color='r')

<Axes: xlabel='x', ylabel='y'>

correlation(nonlinear, 'x', 'y')

np.float64(0.0)

La corrélation est affectée par les valeurs aberrantes

Les valeurs aberrantes peuvent avoir un effet important sur la corrélation. Voici un exemple où un diagramme de dispersion pour lequel $r$ est égal à 1 est transformé en un diagramme pour lequel $r$ est égal à 0, par l'ajout d'un seul point aberrant.

line = pd.DataFrame({
    'x': np.array([1, 2, 3, 4]),
    'y': np.array([1, 2, 3, 4])
})

sns.scatterplot(data=line, x='x', y='y', s=30, color='r')

<Axes: xlabel='x', ylabel='y'>

correlation(line, 'x', 'y')

np.float64(1.0)

outlier = pd.DataFrame({
    'x': np.array([1, 2, 3, 4, 5]),
    'y': np.array([1, 2, 3, 4, 0])
})

sns.scatterplot(data=outlier, x='x', y='y', s=30, color='r')

<Axes: xlabel='x', ylabel='y'>

correlation(outlier, 'x', 'y')

np.float64(0.0)

Les corrélations écologiques doivent être interprétées avec prudence

Les corrélations basées sur des données agrégées peuvent être trompeuses. À titre d'exemple, voici des données sur les résultats du SAT en lecture critique et en mathématiques en 2014. Il y a un point pour chacun des 50 États et un pour Washington, D.C. La colonne "Taux de participation" contient le pourcentage de lycéens de dernière année qui ont passé le test. Les trois colonnes suivantes indiquent le score moyen de l'État pour chaque partie du test, et la dernière colonne est la moyenne des scores totaux du test.

sat2014 = pd.read_csv(path_data + 'sat2014.csv').sort_values('State')
sat2014

	State	Participation Rate	Critical Reading	Math	Writing	Combined
21	Alabama	6.7	547	538	532	1617
34	Alaska	54.2	507	503	475	1485
26	Arizona	36.4	522	525	500	1547
15	Arkansas	4.2	573	571	554	1698
33	California	60.3	498	510	496	1504
12	Colorado	14.3	582	586	567	1735
30	Connecticut	88.4	507	510	508	1525
49	Delaware	100.0	456	459	444	1359
50	District of Columbia	100.0	440	438	431	1309
43	Florida	72.2	491	485	472	1448
44	Georgia	77.2	488	485	472	1445
41	Hawaii	62.6	484	504	472	1460
48	Idaho	100.0	458	456	450	1364
1	Illinois	4.6	599	616	587	1802
38	Indiana	70.5	497	500	477	1474
2	Iowa	3.1	605	611	578	1794
9	Kansas	5.3	591	596	566	1753
10	Kentucky	4.6	589	585	572	1746
18	Louisiana	4.6	561	556	550	1667
47	Maine	95.6	467	471	449	1387
39	Maryland	78.5	492	495	481	1468
24	Massachusetts	84.1	516	531	509	1556
5	Michigan	3.8	593	610	581	1784
4	Minnesota	5.9	598	610	578	1786
13	Mississippi	3.2	583	566	565	1714
7	Missouri	4.2	595	597	579	1771
20	Montana	17.9	555	552	530	1637
11	Nebraska	3.7	589	587	569	1745
42	Nevada	54.2	495	494	469	1458
23	New Hampshire	70.3	524	530	512	1566
29	New Jersey	79.3	501	523	502	1526
22	New Mexico	12.3	548	543	526	1617
40	New York	76.3	488	502	478	1468
35	North Carolina	63.8	499	507	477	1483
0	North Dakota	2.3	612	620	584	1816
19	Ohio	15.1	555	562	535	1652
16	Oklahoma	4.5	576	571	550	1697
27	Oregon	47.9	523	522	499	1544
36	Pennsylvania	71.4	497	504	480	1481
37	Rhode Island	73.0	497	496	487	1480
45	South Carolina	64.9	488	490	465	1443
3	South Dakota	2.9	604	609	579	1792
14	Tennessee	7.9	578	570	566	1714
46	Texas	62.0	476	495	461	1432
17	Utah	5.2	571	568	551	1690
25	Vermont	63.1	522	525	507	1554
28	Virginia	73.1	518	515	497	1530
32	Washington	63.1	510	518	491	1519
31	West Virginia	14.8	517	505	500	1522
6	Wisconsin	3.9	596	608	578	1782
8	Wyoming	3.3	590	599	573	1762

Le diagramme de dispersion des résultats en mathématiques par rapport aux résultats en lecture critique est très étroitement groupé autour d'une ligne droite ; la corrélation est proche de 0,985.

sns.scatterplot(data=sat2014, x='Critical Reading', y='Math')

<Axes: xlabel='Critical Reading', ylabel='Math'>

correlation(sat2014, 'Critical Reading', 'Math')

np.float64(0.9847558411067434)

Il s'agit d'une corrélation extrêmement élevée. Mais il est important de noter que cela ne reflète pas la force de la relation entre les résultats en mathématiques et en lecture critique des élèves.

Les données consistent en des résultats moyens dans chaque État. Or, ce ne sont pas les États qui passent les tests, mais les élèves. Les données du tableau ont été créées en regroupant tous les étudiants de chaque État en un seul point correspondant aux valeurs moyennes des deux variables dans cet État. Mais tous les élèves de l'État ne se situent pas à ce point, car leurs performances varient. Si vous tracez un point pour chaque élève au lieu d'un seul pour chaque État, il y aura un nuage de points autour de chaque point dans la figure ci-dessus. L'image globale sera plus floue. La corrélation entre les résultats en mathématiques et en lecture critique des élèves sera inférieure à la valeur calculée sur la base des moyennes des États.

Les corrélations basées sur des agrégats et des moyennes sont appelées "corrélations écologiques" et sont fréquemment rapportées. Comme nous venons de le voir, elles doivent être interprétées avec prudence.

Sérieux ou pince-sans-rire ?

En 2012, un [article] (http://www.biostat.jhsph.edu/courses/bio621/misc/Chocolate%20consumption%20cognitive%20function%20and%20nobel%20laurates%20%28NEJM%29.pdf) publié dans le très respecté New England Journal of Medicine a examiné la relation entre la consommation de chocolat et les prix Nobel dans un groupe de pays. Le Scientific American a réagi sérieusement alors que Les autres ont été plus détendus. Vous êtes invités à prendre votre propre décision ! Le graphique suivant, fourni dans l'article, devrait vous inciter à aller y jeter un coup d'œil.

from IPython.display import Image
Image("images/chocoNobel.png")

La ligne de régression

Le coefficient de corrélation $r$ ne mesure pas seulement le degré de regroupement des points d'un nuage de points autour d'une ligne droite. Il permet également d'identifier la ligne droite autour de laquelle les points sont regroupés. Dans cette section, nous allons retracer le chemin parcouru par Galton et Pearson pour découvrir cette ligne.

Voici un ensemble de données historiques utilisées pour la prédiction de la taille des adultes sur la base de la taille de leurs parents. Nous avons étudié cet ensemble de données dans une section précédente. La table hauteurs contient des données sur la taille moyenne des parents et la taille de l'enfant (toutes en pouces) pour une population de 934 "enfants" adultes. Rappelons que la taille moyenne des parents est une moyenne des tailles des deux parents.

# Data on heights of parents and their adult children
original = pd.read_csv(path_data + 'family_heights.csv')

heights = pd.DataFrame({
    'MidParent': original['midparentHeight'],
    'Child': original['childHeight']
})

heights

	MidParent	Child
0	75.43	73.2
1	75.43	69.2
2	75.43	69.0
3	75.43	69.0
4	73.66	73.5
...	...	...
929	66.64	64.0
930	66.64	62.0
931	66.64	61.0
932	65.27	66.5
933	65.27	57.0

934 rows × 2 columns

sns.scatterplot(data=heights, x='MidParent', y='Child')

<Axes: xlabel='MidParent', ylabel='Child'>

L'une des principales raisons de l'étude des données était de pouvoir prédire la taille à l'âge adulte d'un enfant né de parents similaires à ceux de l'ensemble de données.

En regardant le graphe, on peut commencer par regarder la forme du nuage du point qui indique une association positive entre les deux variables: lorsqu'une des variables augmente, l'autre à tendance aussi à augmenter.

Comment utiliser ces données pour prédire la taille d'une nouvelle personne étant donné la taille de ses parents ?

Une idée naturelle est de base notre prédiction sur tous les points dans notre jeu de données pour lesquels la taille moyenne des parents est proche de la taille moyenne des parents de la nouvelle personne.

Pour ce faire, nous avons écrit une fonction appelée predict_child, ci-dessous, qui prend la taille moyenne des parents (pour une nouvelle personne) comme argument et renvoie la taille moyenne de tous les enfants dont la taille moyenne des parents se situe à moins d'un demi-pouce de l'argument dans notre jeu de données.

def predict_child(mpht):
    """Return a prediction of the height of a child 
    whose parents have a midparent height of mpht.
    
    The prediction is the average height of the children 
    whose midparent height is in the range mpht plus or minus 0.5 inches.
    """
    
    close_points = heights[(heights['MidParent'] >= mpht - 0.5) & (heights['MidParent'] <= mpht + 0.5)]
    return close_points['Child'].mean()

Nous avons appliqué la fonction à la colonne des hauteurs des "parents moyens" et visualisé le résultat.

heights['Prediction'] = heights['MidParent'].apply(predict_child)
heights_with_predictions = heights

sns.scatterplot(data=heights_with_predictions, x='MidParent', y='Child')
sns.scatterplot(data=heights_with_predictions, x='MidParent', y='Prediction', color='red')

<Axes: xlabel='MidParent', ylabel='Child'>

La prédiction à une hauteur moyenne donnée se situe approximativement au centre de la bande verticale de points à la hauteur donnée.

Voyons si nous pouvons trouver un moyen d'identifier cette ligne par une formule explicite. Tout d'abord, remarquez que l'association linéaire ne dépend pas des unités de mesure - nous pourrions tout aussi bien mesurer les deux variables en unités standard.

def standard_units(xyz):
    "Convert any array of numbers to standard units."
    return (xyz - np.mean(xyz)) / np.std(xyz)

heights_SU = pd.DataFrame({
    'MidParent SU': standard_units(heights['MidParent']),
    'Child SU': standard_units(heights['Child'])
})
heights_SU

	MidParent SU	Child SU
0	3.454652	1.804156
1	3.454652	0.686005
2	3.454652	0.630097
3	3.454652	0.630097
4	2.472085	1.888017
...	...	...
929	-1.424873	-0.767591
930	-1.424873	-1.326667
931	-1.424873	-1.606205
932	-2.185390	-0.068747
933	-2.185390	-2.724356

934 rows × 2 columns

Sur cette échelle, nous pouvons calculer nos prédictions exactement comme auparavant. Mais nous devons d'abord déterminer comment convertir notre ancienne définition des points "proches" en une valeur sur la nouvelle échelle. Nous avions dit que les hauteurs moyennes des parents étaient "proches" si elles se trouvaient à moins de 0,5 pouce l'une de l'autre. Étant donné que les unités standard mesurent les distances en unités d'écart-type, nous devons déterminer combien d'écarts-type de la taille des parents correspondent à 0,5 pouce.

Un écart-type de la taille des parents est d'environ 1,8 pouce. Ainsi, 0,5 pouce correspond à environ 0,28 écart-type.

sd_midparent = np.std(heights['MidParent'])
sd_midparent

np.float64(1.8014050969207571)

0.5/sd_midparent

np.float64(0.277561110965367)

Nous sommes maintenant prêts à modifier notre fonction de prédiction pour faire des prédictions sur l'échelle des unités standard. Tout ce qui a changé, c'est que nous utilisons le tableau des valeurs en unités standard et que nous définissons "proche" comme ci-dessus.

def predict_child_su(mpht_su):
    """Return a prediction of the height (in standard units) of a child 
    whose parents have a midparent height of mpht_su in standard units.
    """
    close = 0.5 / sd_midparent
    close_points = heights_SU[
        (heights_SU['MidParent SU'] >= mpht_su - close) & 
        (heights_SU['MidParent SU'] <= mpht_su + close)
    ]
    return close_points['Child SU'].mean()

heights_SU['Prediction SU'] = heights_SU['MidParent SU'].apply(predict_child_su)
heights_with_su_predictions = heights_SU

sns.scatterplot(data=heights_with_su_predictions, x='MidParent SU', y='Child SU')
sns.scatterplot(data=heights_with_su_predictions, x='MidParent SU', y='Prediction SU', color='red')

<Axes: xlabel='MidParent SU', ylabel='Child SU'>

Ce graphique ressemble exactement au graphique dessiné à l'échelle originale. Seuls les nombres sur les axes ont changé. Cela confirme que nous pouvons comprendre le processus de prédiction en travaillant simplement avec des unités standard.

Identifier la ligne en unités standard

Le diagramme de dispersion ci-dessus a la forme d'un ballon de football, c'est-à-dire qu'il est à peu près ovale comme un ballon de rugby. Tous les diagrammes de dispersion n'ont pas cette forme, même ceux qui montrent une association linéaire. Mais dans cette section, nous ne travaillerons qu'avec des diagrammes de dispersion en forme de ballon de football. Dans la section suivante, nous généraliserons notre analyse à d'autres formes de diagrammes.

Voici un diagramme de dispersion en forme de ballon de football dont les deux variables sont mesurées en unités standard. La ligne à 45 degrés est représentée en rouge.

r = 0.5
x_demo = np.random.normal(0, 1, 10000)
z_demo = np.random.normal(0, 1, 10000)
y_demo = r * x_demo + np.sqrt(1 - r**2) * z_demo

plots.figure(figsize=(6,6))
plots.xlim(-4, 4)
plots.ylim(-4, 4)
plots.scatter(x_demo, y_demo, s=10)
plots.plot([-4, 4], [-4, 4], color='r', lw=2)
plots.xlabel('x in standard units')
plots.ylabel('y in standard units')

Text(0, 0.5, 'y in standard units')

Mais la ligne des 45 degrés n'est pas la ligne qui détermine les centres des bandes verticales. Vous pouvez le constater dans la figure ci-dessous, où la ligne verticale à 1,5 unité standard est représentée en noir. Les points du nuage de points situés près de la ligne noire ont tous une hauteur comprise entre -2 et 3. La ligne rouge est trop haute pour que l'on puisse en déterminer le centre.

r = 0.5
x_demo = np.random.normal(0, 1, 10000)
z_demo = np.random.normal(0, 1, 10000)
y_demo = r * x_demo + np.sqrt(1 - r**2) * z_demo

plots.figure(figsize=(6,6))
plots.xlim(-4, 4)
plots.ylim(-4, 4)
plots.scatter(x_demo, y_demo, s=10)
plots.plot([-4, 4], [-4, 4], color='r', lw=2)
plots.plot([1.5, 1.5], [-4, 4], color='k', lw=2)
plots.xlabel('x in standard units')
plots.ylabel('y in standard units')

Text(0, 0.5, 'y in standard units')

La ligne à 45 degrés n'est donc pas le "graphique des moyennes". Cette ligne est la ligne verte illustrée ci-dessous.

r = 0.5
x_demo = np.random.normal(0, 1, 10000)
z_demo = np.random.normal(0, 1, 10000)
y_demo = r*x_demo + np.sqrt(1 - r**2)*z_demo
plots.figure(figsize=(6,6))
plots.xlim(-4, 4)
plots.ylim(-4, 4)
plots.scatter(x_demo, y_demo, s=10)
plots.plot([-4, 4], [-4*0.6,4*0.6], color='g', lw=2)
plots.plot([-4,4],[-4,4], color='r', lw=2)
plots.plot([1.5,1.5], [-4,4], color='k', lw=2)
plots.xlabel('x in standard units')
plots.ylabel('y in standard units')

Text(0, 0.5, 'y in standard units')

Les deux lignes passent par l'origine (0, 0). La ligne verte passe par le centre des bandes verticales (du moins grossièrement) et est plus plate que la ligne rouge à 45 degrés.

La pente de la ligne à 45 degrés est de 1. La pente de la ligne verte du "graphique des moyennes" est donc une valeur positive mais inférieure à 1.

De quelle valeur s'agit-il ? Vous l'avez deviné, c'est $r$.

La ligne de régression, en unités standard

La ligne verte du "graphique des moyennes" est appelée "ligne de régression", pour une raison que nous expliquerons bientôt. Mais d'abord, simulons quelques diagrammes de dispersion en forme de ballon de football avec différentes valeurs de $r$, et voyons comment la ligne change. Dans chaque cas, la ligne rouge à 45 degrés a été tracée à des fins de comparaison.

La fonction qui effectue la simulation s'appelle regression_line et prend $r$ comme argument.

def regression_line(r):
    x = np.random.normal(0, 1, 10000)
    z = np.random.normal(0, 1, 10000)
    y = r * x + (np.sqrt(1 - r**2)) * z
    plots.figure(figsize=(6, 6))
    plots.xlim(-4, 4)
    plots.ylim(-4, 4)
    plots.scatter(x, y)
    plots.plot([-4, 4], [-4 * r, 4 * r], color='g', lw=2)
    if r >= 0:
        plots.plot([-4, 4], [-4, 4], lw=2, color='r')
    else:
        plots.plot([-4, 4], [4, -4], lw=2, color='r')

regression_line(0.95)

regression_line(0.6)

Lorsque $r$ est proche de 1, le diagramme de dispersion, la ligne à 45 degrés et la ligne de régression sont tous très proches les uns des autres. Mais pour des valeurs plus modérées de $r$, la ligne de régression est nettement plus plate.

L'effet de régression

En termes de prédiction, cela signifie que pour un parent dont la taille médiane est de 1,5 unité standard, notre prédiction de la taille de l'enfant est un peu moins que 1,5 unité standard. Si la taille médiane des parents est de 2 unités standard, nous prédisons que la taille de l'enfant sera légèrement inférieure à 2 unités standard.

En d'autres termes, nous prédisons que l'enfant sera un peu plus proche de la moyenne que ne l'étaient ses parents. C'est ce qu'on appelle la "régression à la moyenne" et c'est ainsi que le nom régression est apparu.

La régression à la moyenne fonctionne également lorsque la taille moyenne des parents est inférieure à la moyenne. Les enfants dont la taille moyenne des parents était inférieure à la moyenne se sont avérés être un peu plus grands par rapport à leur génération, en moyenne.

En général, on s'attend à ce que les individus qui s'éloignent de la moyenne sur une variable ne s'éloignent pas autant de la moyenne sur l'autre variable. C'est ce qu'on appelle l'"effet de régression".

Gardez à l'esprit que l'effet de régression est une affirmation sur les moyennes. Par exemple, il indique que si vous prenez tous les enfants dont la taille moyenne des parents est de 1,5 unité standard, la taille moyenne de ces enfants est légèrement inférieure à 1,5 unité standard. Cela ne signifie pas que tous ces enfants auront une taille légèrement inférieure à 1,5 unité standard. Certains seront plus grands, d'autres plus petits. La moyenne de ces tailles sera inférieure à 1,5 unité standard.

L'équation de la droite de régression

Dans la régression, nous utilisons la valeur d'une variable (que nous appellerons $x$) pour prédire la valeur d'une autre (que nous appellerons $y$). Lorsque les variables $x$ et $y$ sont mesurées en unités standard, la droite de régression permettant de prédire $y$ en fonction de $x$ a une pente de $r$ et passe par l'origine. L'équation de la droite de régression peut donc s'écrire comme suit :

$$ \mbox{estimation de }y ~=~ r \cdot x ~~~ \mbox{lorsque les deux variables sont mesurées en unités standard} $$

Dans les unités originales des données, cela devient

$$ \frac{\mbox{estimation de}~y ~-~\mbox{moyenne de}~y}{\mbox{SD de}~y} ~=~ r \times \frac{\mbox{la donnée}~x ~-~\mbox{moyenne de}~x}{\mbox{SD de}~x} $$

La pente et l'ordonnée à l'origine de la droite de régression peuvent être déduites du diagramme ci-dessus.

$$ \mathbf{\mbox{pente de la droite de régression}} ~=~ r \cdot \frac{\mbox{SD de }y}{\mbox{SD de }x} $$

$$ \mathbf{\mbox{interception de la droite de régression}} ~=~ \mbox{moyenne de }y ~-~ \mbox{pente} \cdot \mbox{moyenne de }x $$

Les trois fonctions ci-dessous calculent la corrélation, la pente et l'ordonnée à l'origine. Elles prennent toutes trois arguments : le nom du tableau, l'intitulé de la colonne contenant $x$ et l'intitulé de la colonne contenant $y$.

def correlation(t, label_x, label_y):
    return np.mean(standard_units(t[label_x]) * standard_units(t[label_y]))

def slope(t, label_x, label_y):
    r = correlation(t, label_x, label_y)
    return r * np.std(t[label_y]) / np.std(t[label_x])

def intercept(t, label_x, label_y):
    return np.mean(t[label_y]) - slope(t, label_x, label_y) * np.mean(t[label_x])

La ligne de régression dans les unités des données

La corrélation entre la taille moyenne des parents et la taille de l'enfant est de 0,32 :

family_r = correlation(heights, 'MidParent', 'Child')
family_r

np.float64(0.32094989606395924)

Nous pouvons également trouver l'équation de la droite de régression pour prédire la taille de l'enfant en fonction de la taille moyenne des parents.

family_slope = slope(heights, 'MidParent', 'Child')
family_intercept = intercept(heights, 'MidParent', 'Child')
family_slope, family_intercept

(np.float64(0.637360896969479), np.float64(22.63624054958975))

L'équation de la droite de régression est la suivante

$$ \mbox{estimation de la taille de l'enfant} ~=~ 0.64 \cdot \mbox{taille du parent moyen} ~+~ 22.64 $$

Cette équation est également connue sous le nom d'équation de régression. La principale utilisation de l'équation de régression est de prédire $y$ sur la base de $x$.

Par exemple, pour une taille moyenne des parents de 70,48 pouces, l'équation de régression prédit que la taille de l'enfant sera de 67,56 pouces.

family_slope * 70.48 + family_intercept

np.float64(67.55743656799862)

Notre prédiction initiale, créée en prenant la taille moyenne de tous les enfants dont la taille moyenne des parents était proche de 70,48, s'est avérée assez proche : 67,63 pouces par rapport à la prédiction de la ligne de régression de 67,55 pouces.

heights_with_predictions[heights_with_predictions['MidParent'] == 70.48].head(3)

	MidParent	Child	Prediction
33	70.48	74.0	67.634239
34	70.48	70.0	67.634239
35	70.48	68.0	67.634239

Voici toutes les lignes du tableau, ainsi que nos prédictions initiales et les nouvelles prédictions de régression de la taille des enfants.

heights_with_predictions['Regression Prediction'] = family_slope * heights_with_predictions['MidParent'] + family_intercept
heights_with_predictions

	MidParent	Child	Prediction	Regression Prediction
0	75.43	73.2	70.100000	70.712373
1	75.43	69.2	70.100000	70.712373
2	75.43	69.0	70.100000	70.712373
3	75.43	69.0	70.100000	70.712373
4	73.66	73.5	70.415789	69.584244
...	...	...	...	...
929	66.64	64.0	65.156579	65.109971
930	66.64	62.0	65.156579	65.109971
931	66.64	61.0	65.156579	65.109971
932	65.27	66.5	64.229630	64.236786
933	65.27	57.0	64.229630	64.236786

934 rows × 4 columns

sns.scatterplot(data=heights_with_predictions, x='MidParent', y='Child')
sns.scatterplot(data=heights_with_predictions, x='MidParent', y='Prediction', color='red')
sns.scatterplot(data=heights_with_predictions, x='MidParent', y='Regression Prediction', color='green')

<Axes: xlabel='MidParent', ylabel='Child'>

Les points verts indiquent les prédictions de régression, toutes sur la ligne de régression. Remarquez que la ligne est très proche du graphique rouges des moyennes. Pour ces données, la ligne de régression fait un bon travail d'approximation des centres des bandes verticales.

Valeurs ajustées

Les prédictions se situent toutes sur la ligne et sont appelées "valeurs ajustées". La fonction fit prend le nom du tableau et les étiquettes $x$ et $y$, et renvoie un tableau de valeurs ajustées, une valeur ajustée pour chaque point du nuage de points.

def fit(table, x, y):
    """Return the height of the regression line at each x value."""
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

Il est plus facile de voir la ligne dans le graphique ci-dessous que dans le graphique ci-dessus.

heights['Fitted'] = fit(heights, 'MidParent', 'Child')
sns.scatterplot(data=heights, x='MidParent', y='Child')
sns.scatterplot(data=heights, x='MidParent', y='Fitted', color='red')

<Axes: xlabel='MidParent', ylabel='Child'>

Une autre façon de tracer la ligne est d'utiliser la fonction lineplot de seaborn (notez que sns.scatterplotet sns.lineplot sont similaires à, respectivement, sns.relplot(..., kind='scatter') et sns.relplot(..., kind='line').

sns.scatterplot(data=heights, x='MidParent', y='Child')
sns.lineplot(data=heights, x='MidParent', y='Fitted', color='red')

<AxesSubplot:xlabel='MidParent', ylabel='Child'>

Unités de mesure de la pente

La pente est un rapport, et il vaut la peine de prendre un moment pour étudier les unités dans lesquelles elle est mesurée. Notre exemple provient d'un ensemble de données concernant les mères ayant accouché dans un système hospitalier. Le nuage de points des poids de grossesse par rapport aux tailles ressemble à un ballon de football qui a été utilisé dans un match de trop, mais il est suffisamment proche d'un ballon de football pour que nous puissions justifier le fait d'y placer notre ligne ajustée. Dans les sections suivantes, nous verrons comment rendre ces justifications plus formelles.

baby = pd.read_csv(path_data + 'baby.csv')

sns.scatterplot(data=baby, x='Maternal Height', y='Maternal Pregnancy Weight')
a = slope(baby, 'Maternal Height', 'Maternal Pregnancy Weight')
b = intercept(baby, 'Maternal Height', 'Maternal Pregnancy Weight')
sns.lineplot(x=baby['Maternal Height'], y=a * baby['Maternal Height'] + b, color='red')

<Axes: xlabel='Maternal Height', ylabel='Maternal Pregnancy Weight'>

slope(baby, 'Maternal Height', 'Maternal Pregnancy Weight')

np.float64(3.572846259275056)

La pente de la ligne de régression est 3,57 livres par pouce. Cela signifie que pour deux femmes dont la taille diffère d'un pouce, notre prédiction du poids de la grossesse sera différente de 3,57 livres. Pour une femme qui mesure 2 pouces de plus qu'une autre (soit 2*2.54=5.08 cm), notre prédiction du poids de grossesse sera de

$$ 2 \times 3,57 ~=~ 7,14 $$

livres (soit $0.453592*7.14 \approx 3.24$kg) de plus que notre prédiction pour la femme la plus petite.

Remarquez que les bandes verticales successives du nuage de points sont espacées d'un pouce, car les hauteurs ont été arrondies au pouce le plus proche. Une autre façon de considérer la pente est de prendre deux bandes consécutives (qui sont nécessairement espacées d'un pouce), correspondant à deux groupes de femmes séparées par un pouce de hauteur. La pente de 3,57 livres par pouce signifie que le poids de grossesse moyen du groupe le plus grand est supérieur d'environ 3,57 livres à celui du groupe le plus petit.

Exemple

Supposons que notre objectif soit d'utiliser la régression pour estimer la taille d'un basset en fonction de son poids, en utilisant un échantillon qui semble compatible avec le modèle de régression. Supposons que la corrélation observée $r$ soit de 0,5 et que les statistiques sommaires pour les deux variables soient celles du tableau ci-dessous :

Variable	Moyenne	Écart-type
Taille	14 pouces	2 pouces
Poids	50 livres	5 livres

Pour calculer l'équation de la droite de régression, nous avons besoin de la pente et de l'ordonnée à l'origine.

$$ \mbox{slope} ~=~ \frac{r \cdot \mbox{SD de }y}{\mbox{SD de }x} ~=~ \frac{0.5 \cdot 2 \mbox{ pouces}}{5 \mbox{ livres}} ~=~ 0.2 ~\mbox{pouces par livre} $$

$$ \mbox{intercept} ~=~ \mbox{moyenne de }y - \mbox{pente}\cdot \mbox{moyenne de } x ~=~ 14 \mbox{ pouces} ~-~ 0.2 \mbox{ pouces par livre} \cdot 50 \mbox{ livres} ~=~ 4 \mbox{ pouces} $$

L'équation de la droite de régression nous permet de calculer la taille estimée, en pouces, à partir d'un poids donné en livres :

$$ \mbox{taille estimée} ~=~ 0.2 \cdot \mbox{poids donné} ~+~ 4 $$

La pente de la ligne mesure l'augmentation de la taille estimée par unité d'augmentation du poids. La pente est positive, et il est important de noter que cela ne signifie pas que nous pensons que les bassets deviennent plus grands s'ils prennent du poids. La pente reflète la différence entre les tailles moyennes de deux groupes de chiens dont le poids est différent d'une livre. Plus précisément, considérons un groupe de chiens dont le poids est de $w$ livres et un groupe dont le poids est de $w+1$ livres. On estime que le second groupe est plus grand de 0,2 pouce, en moyenne. Ceci est vrai pour toutes les valeurs de $w$ dans l'échantillon.

En général, la pente de la ligne de régression peut être interprétée comme l'augmentation moyenne de $y$ par unité d'augmentation de $x$. Notez que si la pente est négative, alors pour chaque augmentation unitaire de $x$, la moyenne de $y$ diminue.

Note de fin

Même si nous n'établirons pas la base mathématique de l'équation de régression, nous pouvons voir qu'elle donne d'assez bonnes prédictions lorsque le nuage de points a la forme d'un ballon de football. C'est un fait mathématique surprenant que, quelle que soit la forme du nuage de points, la même équation donne la "meilleure" parmi toutes les lignes droites. C'est le sujet de la section suivante.

Retour à l'exemple du roman "Little women"

Nous avons dit précédemment que la formule que nous venons d'obtenir donne exactement le même résultat que la méthode de régression aux moindres carrés vue précédemment. Cela n'a rien d'évident a priori. Nous n'allons pas démontrer mathématiquement que c'est vrai, mais nous pouvons au moins vérifier empiriquement que c'est le cas sur l'exemple du roman "Little women".

def standard_units(any_numbers):
    "Convert any array of numbers to standard units."
    return (any_numbers - np.mean(any_numbers)) / np.std(any_numbers)

def correlation(t, x, y):
    return np.mean(standard_units(t[x]) * standard_units(t[y]))

def slope(table, x, y):
    r = correlation(table, x, y)
    return r * np.std(table[y]) / np.std(table[x])

def intercept(table, x, y):
    a = slope(table, x, y)
    return np.mean(table[y]) - a * np.mean(table[x])

def fit(table, x, y):
    """Return the height of the regression line at each x value."""
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

correlation(little_women, 'Periods', 'Characters')

np.float64(0.9229576895854816)

Le diagramme de dispersion est remarquablement proche de la linéarité et la corrélation est supérieure à 0,92.

little_women['Linear Prediction'] = fit(little_women, 'Periods', 'Characters')
sns.scatterplot(data=little_women, x='Periods', y='Characters')
sns.lineplot(data=little_women, x='Periods', y='Linear Prediction', color='red')

<Axes: xlabel='Periods', ylabel='Characters'>

actual = little_women['Characters']
predicted = little_women['Linear Prediction']
errors = actual - predicted

little_women['Error'] = errors
lw_with_predictions = little_women
lw_with_predictions

	Periods	Characters	Linear Prediction	Error
0	189	21759	21183.596794	575.403206
1	188	22148	21096.618953	1051.381047
2	231	20558	24836.666127	-4278.666127
3	195	25526	21705.463842	3820.536158
4	255	23395	26924.134317	-3529.134317
5	140	14622	16921.682573	-2299.682573
6	131	14431	16138.882001	-1707.882001
7	214	22476	23358.042826	-882.042826
8	337	33767	34056.317301	-289.317301
9	185	18508	20835.685429	-2327.685429
10	193	23024	21531.508159	1492.491841
11	429	39363	42058.278696	-2695.278696
12	175	18669	19965.907017	-1296.907017
13	180	17924	20400.796223	-2476.796223
14	181	17385	20487.774064	-3102.774064
15	172	15353	19704.973493	-4351.973493
16	155	13879	18226.350192	-4347.350192
17	144	17432	17269.593938	162.406062
18	121	17338	15269.103589	2068.896411
19	145	14252	17356.571779	-3104.571779
20	269	23545	28141.824095	-4596.824095
21	120	13708	15182.125748	-1474.125748
22	247	23015	26228.311587	-3213.311587
23	182	24239	20574.751906	3664.248094
24	91	13526	12659.768351	866.231649
25	150	21739	17791.460985	3947.539015
26	109	15754	14225.369494	1528.630506
27	271	30992	28315.779778	2676.220222
28	233	26254	25010.621810	1243.378190
29	218	24137	23705.954191	431.045809
30	178	20723	20226.840541	496.159459
31	224	23158	24227.821238	-1069.821238
32	232	25777	24923.643969	853.356031
33	257	32496	27098.090000	5397.910000
34	201	21310	22227.330889	-917.330889
35	100	11441	13442.568922	-2001.568922
36	157	23524	18400.305874	5123.694126
37	206	26091	22662.220096	3428.779904
38	263	27473	27619.957048	-146.957048
39	61	11368	10050.433113	1317.566887
40	187	26464	21009.641112	5454.358888
41	118	16753	15008.170065	1744.829935
42	305	33202	31273.026380	1928.973620
43	95	10289	13007.679716	-2718.679716
44	96	12558	13094.657557	-536.657557
45	234	27094	25097.599651	1996.400349
46	392	40935	38840.098570	2094.901430

Nous pouvons utiliser slope et intercept pour calculer la pente et l'ordonnée à l'origine de la droite ajustée.

lw_reg_slope = slope(little_women, 'Periods', 'Characters')
lw_reg_intercept = intercept(little_women, 'Periods', 'Characters')

print('Slope of Regression Line:    ', np.round(lw_reg_slope), 'characters per period')
print('Intercept of Regression Line:', np.round(lw_reg_intercept), 'characters')
lw_errors(lw_reg_slope, lw_reg_intercept)

Slope of Regression Line:     87.0 characters per period
Intercept of Regression Line: 4745.0 characters

On retrouve bien les valeurs obtenues par la méthodes des moindres carrés !

Régression par les moindres carrés, plus d'exemples (facultatif)

Dans une section précédente, nous avons développé des formules pour la pente et l'ordonnée à l'origine de la droite de régression à l'aide d'un diagramme de dispersion en forme de ballon de football. Il s'avère que la pente et l'ordonnée à l'origine de la droite des moindres carrés ont les mêmes formules que celles que nous avons développées, quelle que soit la forme du diagramme de dispersion.

Nous l'avons vu dans l'exemple de Little Women, mais confirmons-le dans un exemple où le diagramme de dispersion n'a manifestement pas la forme d'un ballon de football. Pour les données, nous sommes une fois de plus redevables aux riches [archives de données du professeur Larry Winner] (http://www.stat.ufl.edu/~winner/datasets.html) de l'université de Floride. Une étude de 2013 parue dans l'International Journal of Exercise Science a étudié des athlètes collégiaux de lancer de poids et a examiné la relation entre la force et la distance du lancer de poids. La population est composée de 28 athlètes féminines de niveau collégial. La force a été mesurée par la plus grande quantité (en kilogrammes) que l'athlète a soulevée dans le "1RM power clean" pendant la pré-saison. La distance (en mètres) correspondait au record personnel de l'athlète.

def standard_units(any_numbers):
    "Convert any array of numbers to standard units."
    return (any_numbers - np.mean(any_numbers)) / np.std(any_numbers)

def correlation(t, x, y):
    return np.mean(standard_units(t[x]) * standard_units(t[y]))

def slope(table, x, y):
    r = correlation(table, x, y)
    return r * np.std(table[y]) / np.std(table[x])

def intercept(table, x, y):
    a = slope(table, x, y)
    return np.mean(table[y]) - a * np.mean(table[x])

def fit(table, x, y):
    """Return the height of the regression line at each x value."""
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

shotput = pd.read_csv(path_data + 'shotput.csv')

shotput

	Weight Lifted	Shot Put Distance
0	37.5	6.4
1	51.5	10.2
2	61.3	12.4
3	61.3	13.0
4	63.6	13.2
5	66.1	13.0
6	70.0	12.7
7	92.7	13.9
8	90.5	15.5
9	90.5	15.8
10	94.8	15.8
11	97.0	16.8
12	97.0	17.1
13	97.0	17.8
14	102.0	14.8
15	102.0	15.5
16	103.6	16.1
17	100.4	16.2
18	108.4	17.9
19	114.0	15.9
20	115.3	15.8
21	114.9	16.7
22	114.7	17.6
23	123.6	16.8
24	125.8	17.0
25	119.1	18.2
26	118.9	19.2
27	141.1	18.6

sns.scatterplot(data=shotput, x='Weight Lifted', y='Shot Put Distance')

<AxesSubplot:xlabel='Weight Lifted', ylabel='Shot Put Distance'>

Il ne s'agit pas d'un diagramme de dispersion en forme de ballon de football. En fait, il semble avoir une légère composante non linéaire. Mais si nous insistons pour utiliser une ligne droite pour faire nos prédictions, il y a toujours une meilleure ligne droite parmi toutes les lignes droites.

Nos formules pour la pente et l'ordonnée à l'origine de la droite de régression, dérivées des diagrammes de dispersion en forme de ballon de football, donnent les valeurs suivantes.

slope(shotput, 'Weight Lifted', 'Shot Put Distance')

0.09834382159781994

intercept(shotput, 'Weight Lifted', 'Shot Put Distance')

5.959629098373956

Est-il toujours utile d'utiliser ces formules même si le nuage de points n'a pas la forme d'un ballon de football ? Nous pouvons répondre à cette question en trouvant la pente et l'ordonnée à l'origine de la droite qui minimise le mse.

Nous allons définir la fonction shotput_linear_mse pour prendre une pente et une ordonnée à l'origine comme arguments et retourner le mse correspondant. Ensuite, minimize appliqué à shotput_linear_mse renverra la meilleure pente et la meilleure ordonnée à l'origine.

def shotput_linear_mse(any_slope, any_intercept):
    x = shotput['Weight Lifted']
    y = shotput['Shot Put Distance']
    fitted = any_slope * x + any_intercept
    return np.mean((y - fitted) ** 2)

def shotput_mse_for_minimize(params):
    slope, intercept = params
    return shotput_linear_mse(slope, intercept)

result = minimize(shotput_mse_for_minimize, x0=[0, 0])  # Initial guess for slope and intercept
best_shotput_params = result.x  # Extracting the best slope and intercept
best_shotput_params

array([0.09834368, 5.95964276])

Ces valeurs sont les mêmes que celles que nous avons obtenues en utilisant nos formules. En résumé :

Quelle que soit la forme du nuage de points, il existe une ligne unique qui minimise l'erreur quadratique moyenne d'estimation. Elle s'appelle la droite de régression, et sa pente et son ordonnée à l'origine sont données par

$$ \mathbf{\mbox{pente de la droite de régression}} ~=~ r \cdot \frac{\mbox{SD de }y}{\mbox{SD de }x} $$

$$ \mathbf{\mbox{intercept de la droite de régression}} ~=~ \mbox{moyenne de }y ~-~ \mbox{pente} \cdot \mbox{moyenne de }x $$

shotput['Best Straight Line'] = fit(shotput, 'Weight Lifted', 'Shot Put Distance')
sns.scatterplot(data=shotput, x='Weight Lifted', y='Shot Put Distance')
sns.lineplot(data=shotput, x='Weight Lifted', y='Best Straight Line', color='red')

<AxesSubplot:xlabel='Weight Lifted', ylabel='Shot Put Distance'>

Régression non linéaire

Le graphique ci-dessus confirme notre observation précédente, à savoir que le nuage de points est légèrement incurvé. Il est donc préférable d'ajuster une courbe plutôt qu'une ligne droite. L'[étude] (http://digitalcommons.wku.edu/ijes/vol6/iss2/10/) postule une relation quadratique entre le poids soulevé et la distance du lancer de poids. Utilisons donc des fonctions quadratiques comme prédicteurs et voyons si nous pouvons trouver la meilleure.

Nous devons trouver la meilleure fonction quadratique parmi toutes les fonctions quadratiques, plutôt que la meilleure ligne droite parmi toutes les lignes droites. La méthode des moindres carrés nous permet de le faire.

Les mathématiques de cette minimisation sont compliquées et il n'est pas facile de s'en rendre compte en examinant le nuage de points. Mais la minimisation numérique est tout aussi facile qu'avec les variables prédicteurs linéaires ! Nous pouvons obtenir le meilleur prédicteur quadratique en utilisant à nouveau minimize. Voyons comment cela fonctionne.

Rappelons qu'une fonction quadratique a la forme suivante

$$ f(x) ~=~ ax^2 + bx + c $$

pour les constantes $a$, $b$ et $c$.

Pour trouver la meilleure fonction quadratique permettant de prédire la distance en fonction du poids soulevé, en utilisant le critère des moindres carrés, nous allons d'abord écrire une fonction qui prend les trois constantes comme arguments, calcule les valeurs ajustées à l'aide de la fonction quadratique ci-dessus, puis renvoie l'erreur quadratique moyenne.

Cette fonction s'appelle shotput_quadratic_mse. Notez que la définition est analogue à celle de lw_mse, sauf que les valeurs ajustées sont basées sur une fonction quadratique au lieu d'une fonction linéaire.

def shotput_quadratic_mse(a, b, c):
    x = shotput['Weight Lifted']
    y = shotput['Shot Put Distance']
    fitted = a * (x**2) + b * x + c
    return np.mean((y - fitted) ** 2)

Nous pouvons maintenant utiliser minimize comme précédemment pour trouver les constantes qui minimisent l'erreur quadratique moyenne.

def shotput_quadratic_mse_for_minimize(params):
    a, b, c = params
    return shotput_quadratic_mse(a, b, c)

result = minimize(shotput_quadratic_mse_for_minimize, x0=[0, 0, 0])  # Initial guess for a, b, c
best = result.x  # Extracting the best parameters
best

array([-1.04178135e-03,  2.83015120e-01, -1.54429228e+00])

Notre prédiction de la distance du lancer de poids pour un athlète qui soulève $x$ kilogrammes est d'environ

$$ -0,00104x^2 ~+~ 0,2827x - 1,5318 $$

mètres. Par exemple, si l'athlète peut soulever 100 kilogrammes, la distance prédite est de 16,33 mètres. Sur le diagramme de dispersion, cette distance est proche du centre d'une bande verticale autour de 100 kilogrammes.

(-0.00104)*(100**2) + 0.2827*100 - 1.5318

16.3382

Voici les prévisions pour toutes les valeurs de Weight Lifted. Vous pouvez voir qu'elles passent par le centre du nuage de points, de façon approximative.

x = shotput['Weight Lifted']
shotput_fit = best[0] * (x**2) + best[1] * x + best[2]

shotput['Best Quadratic Curve'] = shotput_fit
sns.scatterplot(data=shotput, x='Weight Lifted', y='Shot Put Distance')
sns.lineplot(data=shotput, x='Weight Lifted', y='Best Quadratic Curve', color='red')

<AxesSubplot:xlabel='Weight Lifted', ylabel='Shot Put Distance'>

**Note : **Nous avons ajusté une quadratique ici parce que cela a été suggéré dans l'étude originale. Mais il convient de noter qu'à l'extrémité droite du graphique, la courbe quadratique semble proche d'un pic, après quoi la courbe commencera à descendre. Nous ne devrions donc pas utiliser ce modèle pour les nouveaux athlètes qui peuvent soulever des poids beaucoup plus élevés que ceux de notre ensemble de données.

Diagnostics visuels (facultatif)

Supposons qu'un scientifique des données ait décidé d'utiliser la régression linéaire pour estimer les valeurs d'une variable (appelée variable réponse) en fonction d'une autre variable (appelée prédicteur). Pour connaître les performances de cette méthode d'estimation, le scientifique des données doit mesurer l'écart entre les estimations et les valeurs réelles. Ces différences sont appelées résidus.

$$ \mbox{résidu} ~=~ \mbox{valeur observée} ~-~ \mbox{estimation de régression} $$

Un résidu est ce qui reste - le résidu - après l'estimation.

Les résidus sont les distances verticales des points par rapport à la droite de régression. Il y a un résidu pour chaque point du nuage de points. Le résidu est la différence entre la valeur observée de $y$ et la valeur ajustée de $y$, donc pour le point $(x, y)$,

$$ \mbox{résidu} ~~ = ~~ y ~-~ \mbox{valeur ajustée de }y ~~ = ~~ y ~-~ \mbox{hauteur de la droite de régression à }x $$

La fonction residual calcule les résidus. Le calcul prend en compte toutes les fonctions pertinentes que nous avons déjà définies : unités_standard, corrélation, pente, intercept et fit.

family_heights = pd.read_csv(path_data + 'family_heights.csv')
heights = family_heights[['midparentHeight', 'childHeight']].rename(columns={
    'midparentHeight': 'MidParent',
    'childHeight': 'Child'
})

hybrid = pd.read_csv(path_data + 'hybrid.csv')

def standard_units(x):
    return (x - np.mean(x)) / np.std(x)

def correlation(table, x, y):
    x_in_standard_units = standard_units(table[x])
    y_in_standard_units = standard_units(table[y])
    return np.mean(x_in_standard_units * y_in_standard_units)

def slope(table, x, y):
    r = correlation(table, x, y)
    return r * np.std(table[y]) / np.std(table[x])

def intercept(table, x, y):
    a = slope(table, x, y)
    return np.mean(table[y]) - a * np.mean(table[x])

def fit(table, x, y):
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

def residual(table, x, y):
    return table[y] - fit(table, x, y)

Poursuivant notre exemple d'estimation de la taille des enfants adultes (la réponse) en fonction de la taille moyenne des parents (le prédicteur), calculons les valeurs ajustées et les résidus.

heights['Fitted Value'] = fit(heights, 'MidParent', 'Child')
heights['Residual'] = residual(heights, 'MidParent', 'Child')
heights

	MidParent	Child	Fitted Value	Residual
0	75.43	73.2	70.712373	2.487627
1	75.43	69.2	70.712373	-1.512373
2	75.43	69.0	70.712373	-1.712373
3	75.43	69.0	70.712373	-1.712373
4	73.66	73.5	69.584244	3.915756
...	...	...	...	...
929	66.64	64.0	65.109971	-1.109971
930	66.64	62.0	65.109971	-3.109971
931	66.64	61.0	65.109971	-4.109971
932	65.27	66.5	64.236786	2.263214
933	65.27	57.0	64.236786	-7.236786

934 rows × 4 columns

Lorsqu'il y a tant de variables à traiter, il est toujours utile de commencer par la visualisation. La fonction scatter_fit dessine le diagramme de dispersion des données, ainsi que la ligne de régression.

def scatter_fit(table, x, y):
    sns.scatterplot(data=table, x=x, y=y, s=15)
    sns.lineplot(x=table[x], y=fit(table, x, y), lw=4, color='gold')
    plots.xlabel(x)
    plots.ylabel(y)

scatter_fit(heights, 'MidParent', 'Child')

Un tracé résiduel peut être dessiné en traçant les résidus en fonction de la variable prédictive. C'est ce que fait la fonction residual_plot.

def residual_plot(table, x, y):
    x_array = table[x]
    residuals = residual(table, x, y)
    residuals_table = pd.DataFrame({
        x: x_array,
        'residuals': residuals
    })
    sns.scatterplot(data=residuals_table, x=x, y='residuals', color='red')
    xlims = [min(x_array), max(x_array)]
    plots.plot(xlims, [0, 0], color='darkblue', lw=4)
    plots.title('Residual Plot')

residual_plot(heights, 'MidParent', 'Child')

Les hauteurs des parents moyens sont sur l'axe horizontal, comme dans le diagramme de dispersion original. Mais maintenant, l'axe vertical montre les résidus. Remarquez que le graphique semble être centré autour de la ligne horizontale au niveau 0 (en bleu foncé). Remarquez également que le graphique ne présente aucune tendance à la hausse ou à la baisse. Nous observerons plus tard que cette absence de tendance est vraie pour toutes les régressions.

Diagnostics de régression

Les graphiques des résidus nous aident à évaluer visuellement la qualité d'une analyse de régression linéaire. Ces évaluations sont appelées diagnostics. La fonction regression_diagnostic_plots dessine le diagramme de dispersion original ainsi que le diagramme résiduel pour faciliter la comparaison.

def regression_diagnostic_plots(table, x, y):
    plots.figure()  
    scatter_fit(table, x, y)
    
    plots.figure()  
    residual_plot(table, x, y)

regression_diagnostic_plots(heights, 'MidParent', 'Child')

Ce diagramme des résidus indique que la régression linéaire était une méthode d'estimation raisonnable. Remarquez que les résidus sont distribués de manière assez symétrique au-dessus et au-dessous de la ligne horizontale à 0, ce qui correspond au diagramme de dispersion original qui est à peu près symétrique au-dessus et au-dessous. Remarquez également que la répartition verticale du diagramme est assez uniforme pour les valeurs les plus courantes de la taille des enfants. En d'autres termes, à l'exception de quelques points aberrants, le diagramme n'est pas plus étroit à certains endroits et plus large à d'autres.

En d'autres termes, la précision de la régression semble être à peu près la même sur toute la plage observée de la variable prédictive.

Le tracé des résidus d'une bonne régression ne présente aucune tendance. Les résidus sont à peu près les mêmes, au-dessus et au-dessous de la ligne horizontale à 0, sur toute la plage de la variable prédictive.

Détection de la non-linéarité

Le tracé du diagramme de dispersion des données permet généralement de savoir si la relation entre les deux variables n'est pas linéaire. Toutefois, il est souvent plus facile de repérer la non-linéarité dans un diagramme résiduel que dans le diagramme de dispersion original. Cela s'explique généralement par les échelles des deux diagrammes : le diagramme résiduel nous permet de zoomer sur les erreurs et, par conséquent, de repérer plus facilement les schémas.

Nos données sont un dataset sur l'âge et la longueur des dugongs, qui sont des mammifères marins apparentés aux lamantins et aux vaches de mer (image de Wikimedia Commons). Les données se trouvent dans un tableau appelé dugong. L'âge est mesuré en années et la longueur en mètres. Comme les dugongs n'ont pas l'habitude de noter leur date de naissance, leur âge est estimé sur la base de variables telles que l'état de leurs dents.

dugong = pd.read_csv(path_data + 'dugongs.csv')
dugong = dugong[['Length'] + [col for col in dugong.columns if col != 'Length']]
dugong

	Length	Age
0	1.80	1.0
1	1.85	1.5
2	1.87	1.5
3	1.77	1.5
4	2.02	2.5
5	2.27	4.0
6	2.15	5.0
7	2.26	5.0
8	2.35	7.0
9	2.47	8.0
10	2.19	8.5
11	2.26	9.0
12	2.40	9.5
13	2.39	9.5
14	2.41	10.0
15	2.50	12.0
16	2.32	12.0
17	2.43	13.0
18	2.47	13.0
19	2.56	14.5
20	2.65	15.5
21	2.47	15.5
22	2.64	16.5
23	2.56	17.0
24	2.70	22.5
25	2.72	29.0
26	2.57	31.5

Si nous pouvions mesurer la longueur d'un dugong, que pourrions-nous dire de son âge ? Examinons ce que disent nos données. Voici une régression de l'âge (la réponse) sur la longueur (le prédicteur). La corrélation entre les deux variables est importante (0,83).

correlation(dugong, 'Length', 'Age')

0.8296474554905715

Malgré la forte corrélation, le graphique présente une courbe qui est beaucoup plus visible dans le graphique résiduel.

regression_diagnostic_plots(dugong, 'Length', 'Age')

Si la non-linéarité est perceptible dans le diagramme de dispersion original, elle apparaît plus clairement dans le diagramme résiduel.

À l'extrémité inférieure des longueurs, les résidus sont presque tous positifs ; ensuite, ils sont presque tous négatifs ; puis à nouveau positifs à l'extrémité supérieure des longueurs. En d'autres termes, les estimations de la régression sont trop élevées, puis trop faibles, puis trop élevées. Cela signifie qu'il aurait été préférable d'utiliser une courbe plutôt qu'une ligne droite pour estimer les âges.

**Lorsqu'un tracé résiduel présente un schéma, il se peut qu'il y ait une relation non linéaire entre les variables.

Détection de l'hétéroscédasticité

*L'hétéroscédasticité est un mot qui intéressera certainement ceux qui se préparent aux concours d'orthographe. Pour les scientifiques des données, son intérêt réside dans sa signification, à savoir une "répartition inégale".

Rappelons le tableau "hybride" qui contient des données sur les voitures hybrides aux États-Unis. Voici une régression de l'efficacité énergétique sur le taux d'accélération. L'association est négative : les voitures qui accélèrent rapidement ont tendance à être moins efficaces.

regression_diagnostic_plots(hybrid, 'acceleration', 'mpg')

Remarquez que le tracé des résidus s'évase vers l'extrémité inférieure des accélérations. En d'autres termes, la variabilité de la taille des erreurs est plus importante pour les faibles valeurs d'accélération que pour les valeurs élevées. Une variation inégale est souvent plus facile à remarquer dans un diagramme résiduel que dans le diagramme de dispersion original.

**Si le diagramme résiduel présente une variation inégale autour de la ligne horizontale à 0, les estimations de la régression ne sont pas également précises sur toute la plage de la variable prédictive.

Diagnostics numériques (facultatif)

Outre la visualisation, nous pouvons utiliser les propriétés numériques des résidus pour évaluer la qualité de la régression. Nous ne prouverons pas ces propriétés mathématiquement. Nous allons plutôt les observer par le calcul et voir ce qu'elles nous apprennent sur la régression.

Tous les faits énumérés ci-dessous sont valables pour toutes les formes de diagrammes de dispersion, qu'ils soient linéaires ou non.

family_heights = pd.read_csv(path_data + 'family_heights.csv')
heights = family_heights[['midparentHeight', 'childHeight']].rename(columns={
    'midparentHeight': 'MidParent',
    'childHeight': 'Child'
})

dugong = pd.read_csv(path_data + 'dugongs.csv')
dugong = dugong[['Length'] + [col for col in dugong.columns if col != 'Length']]

hybrid = pd.read_csv(path_data + 'hybrid.csv')

def standard_units(x):
    return (x - np.mean(x)) / np.std(x)

def correlation(table, x, y):
    x_in_standard_units = standard_units(table[x])
    y_in_standard_units = standard_units(table[y])
    return np.mean(x_in_standard_units * y_in_standard_units)

def slope(table, x, y):
    r = correlation(table, x, y)
    return r * np.std(table[y]) / np.std(table[x])

def intercept(table, x, y):
    a = slope(table, x, y)
    return np.mean(table[y]) - a * np.mean(table[x])

def fit(table, x, y):
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

def residual(table, x, y):
    return table[y] - fit(table, x, y)

def scatter_fit(table, x, y):
    sns.scatterplot(data=table, x=x, y=y, s=15)
    sns.lineplot(x=table[x], y=fit(table, x, y), lw=4, color='gold')
    plots.xlabel(x)
    plots.ylabel(y)
    
def residual_plot(table, x, y):
    residuals_table = pd.DataFrame({
        x: table[x],
        'residuals': residual(table, x, y)
    })
    sns.scatterplot(data=residuals_table, x=x, y='residuals', color='red')
    xlims = [min(table[x]), max(table[x])]
    plots.plot(xlims, [0, 0], color='darkblue', lw=4)
    plots.title('Residual Plot')

def regression_diagnostic_plots(table, x, y):
    plots.figure()
    scatter_fit(table, x, y)
    plots.figure()
    residual_plot(table, x, y)

heights['Fitted Value'] = fit(heights, 'MidParent', 'Child')
heights['Residual'] = residual(heights, 'MidParent', 'Child')

Les tracés résiduels ne montrent aucune tendance

Pour chaque régression linéaire, qu'elle soit bonne ou mauvaise, le graphique des résidus ne montre aucune tendance. Dans l'ensemble, il est plat. En d'autres termes, les résidus et la variable prédictive ne sont pas corrélés.

C'est ce que l'on peut voir dans tous les diagrammes de résidus ci-dessus. Nous pouvons également calculer la corrélation entre la variable prédictive et les résidus dans chaque cas.

correlation(heights, 'MidParent', 'Residual')

3.090556599598937e-16

Cela ne ressemble pas à zéro, mais il s'agit d'un tout petit nombre qui est égal à 0 en dehors de l'erreur d'arrondi due au calcul. Le voici à nouveau, avec 10 décimales. Le signe moins est dû à l'erreur d'arrondi ci-dessus.

round(correlation(heights, 'MidParent', 'Residual'), 10)

0.0

dugong['Fitted Value'] = fit(dugong, 'Length', 'Age')
dugong['Residual'] = residual(dugong, 'Length', 'Age')

round(correlation(dugong, 'Length', 'Residual'), 10)

0.0

Moyenne des résidus

Quelle que soit la forme du diagramme de dispersion, la moyenne des résidus est de 0,.

Ceci est analogue au fait que si vous prenez n'importe quelle liste de nombres et que vous calculez la liste des écarts par rapport à la moyenne, la moyenne des écarts est égale à 0.

Dans tous les diagrammes de résidus ci-dessus, vous avez vu la ligne horizontale à 0 passer par le centre du diagramme. Il s'agit d'une visualisation de ce fait.

En guise d'exemple numérique, voici la moyenne des résidus dans la régression des tailles des enfants sur la base des hauteurs moyennes des parents.

round(np.mean(heights['Residual']), 10)

0.0

Il en est de même pour la moyenne des résidus de la régression de l'âge des dugongs sur leur longueur. La moyenne des résidus est égale à 0, sauf erreur d'arrondi.

round(np.mean(dugong['Residual']), 10)

-0.0

L'écart-type des résidus

Quelle que soit la forme du nuage de points, l'écart-type des résidus est une fraction de l'écart-type de la variable réponse. Cette fraction est $\sqrt{1-r^2}$.

$$ \mbox{SD des résidus} ~=~ \sqrt{1 - r^2} \cdot \mbox{SD de }y $$

Nous verrons bientôt comment cela permet de mesurer la précision de l'estimation de la régression. Mais d'abord, confirmons-le par un exemple.

Dans le cas de la taille des enfants et de la taille moyenne des parents, l'écart-type des résidus est d'environ 3,39 pouces.

np.std(heights['Residual'])

3.388079916395343

C'est la même chose que $\sqrt{1-r^2}$ fois l'écart-type de la variable réponse :

r = correlation(heights, 'MidParent', 'Child')
np.sqrt(1 - r**2) * np.std(heights['Child'])

3.388079916395342

Il en va de même pour la régression du kilométrage sur l'accélération des voitures hybrides. La corrélation $r$ est négative (environ -0,5), mais $r^2$ est positif et donc $\sqrt{1-r^2}$ est une fraction.

r = correlation(hybrid, 'acceleration', 'mpg')
r

-0.5060703843771185

r = correlation(hybrid, 'acceleration', 'mpg')

hybrid['fitted mpg'] = fit(hybrid, 'acceleration', 'mpg')
hybrid['residual'] = residual(hybrid, 'acceleration', 'mpg')

np.std(hybrid['residual']), np.sqrt(1 - r**2) * np.std(hybrid['mpg'])

(9.43273683343029, 9.432736833430296)

Voyons maintenant comment l'écart-type des résidus permet de mesurer la qualité de la régression. Rappelons que la moyenne des résidus est égale à 0. Par conséquent, plus l'écart-type des résidus est petit, plus les résidus sont proches de 0. En d'autres termes, si l'écart-type des résidus est petit, la taille globale des erreurs dans la régression est petite.

Les cas extrêmes sont lorsque $r=1$ ou $r=-1$. Dans les deux cas, $\sqrt{1-r^2} = 0$. Par conséquent, les résidus ont une moyenne de 0 et un écart-type de 0 également, et les résidus sont donc tous égaux à 0. La droite de régression fait un travail d'estimation parfait. Comme nous l'avons vu plus tôt dans ce chapitre, si $r = \pm 1$, le nuage de points est une droite parfaite et est identique à la droite de régression, il n'y a donc pas d'erreur dans l'estimation de la régression.

Mais en général, $r$ ne se situe pas aux extrêmes. Si $r$ n'est ni $\pm 1$ ni 0, alors $\sqrt{1-r^2}$ est une fraction correcte, et la taille globale approximative de l'erreur de l'estimation de la régression se situe quelque part entre 0 et la DS de $y$.

Le pire cas est celui où $r = 0$. Dans ce cas, $\sqrt{1-r^2} =1$, et l'écart-type des résidus est égal à l'écart-type de $y$. Ceci est cohérent avec l'observation que si $r=0$, la ligne de régression est une ligne plate à la moyenne de $y$. Dans cette situation, l'erreur quadratique moyenne de la régression est l'écart quadratique moyen par rapport à la moyenne de $y$, qui est l'écart-type de $y$. En pratique, si $r = 0$, il n'y a pas d'association linéaire entre les deux variables, et il n'y a donc aucun avantage à utiliser la régression linéaire.

Une autre façon d'interpréter $r$

Nous pouvons réécrire le résultat ci-dessus pour dire que quelle que soit la forme du nuage de points,

$$ \frac{\mbox{SD des résidus}}{\mbox{SD de }y} ~=~ \sqrt{1-r^2} $$

Un résultat complémentaire est que, quelle que soit la forme du nuage de points, l'écart-type des valeurs ajustées est une fraction de l'écart-type des valeurs observées de $y$. Cette fraction est $\vert r \vert$.

$$ \frac{\mbox{SD des valeurs ajustées}}{\mbox{SD de }y} ~=~ \vert r \vert$. $$

Pour voir où la fraction intervient, il faut remarquer que les valeurs ajustées sont toutes sur la droite de régression alors que les valeurs observées de $y$ sont les hauteurs de tous les points du nuage de points et sont plus variables.

scatter_fit(heights, 'MidParent', 'Child')

Les valeurs ajustées vont de 64 à 71 environ, alors que les tailles de tous les enfants sont beaucoup plus variables, allant de 55 à 80 environ.

Pour vérifier le résultat numériquement, il suffit de calculer les deux côtés de l'identité.

correlation(heights, 'MidParent', 'Child')

0.320949896063959

Il s'agit du rapport entre l'écart-type des valeurs ajustées et l'écart-type des valeurs observées du poids de naissance :

np.std(heights['Fitted Value']) / np.std(heights['Child'])

0.32094989606395896

Le rapport est égal à $r$, ce qui confirme notre résultat.

Où intervient la valeur absolue ? Tout d'abord, il faut noter que les écarts types ne peuvent pas être négatifs, de même qu'un rapport d'écarts types ne peut pas l'être. Que se passe-t-il donc lorsque $r$ est négatif ? L'exemple de l'efficacité énergétique et de l'accélération va nous le montrer.

correlation(hybrid, 'acceleration', 'mpg')

-0.5060703843771185

np.std(hybrid['fitted mpg']) / np.std(hybrid['mpg'])

0.5060703843771186

Le rapport des deux SD est $\vert r \vert$.

Une façon plus standard d'exprimer ce résultat est de rappeler que

$$ \mbox{variance} ~=~ \mbox{écart quadratique moyen par rapport à la moyenne} ~=~ \mbox{SD}^2 $$

et donc en élevant au carré les deux côtés de notre résultat,

$$ \frac{\mbox{variance des valeurs ajustées}}{\mbox{variance de }y} ~=~ r^2 $$

2. Inférence pour la régression (facultatif)

Jusqu'à présent, notre analyse de la relation entre les variables a été purement descriptive. Nous savons comment trouver la meilleure ligne droite à tracer dans un nuage de points. Cette droite est la meilleure dans la mesure où elle présente l'erreur quadratique moyenne d'estimation la plus faible parmi toutes les droites.

Mais que se passerait-il si nos données n'étaient qu'un échantillon d'une population plus large ? Si, dans l'échantillon, nous avons trouvé une relation linéaire entre les deux variables, en serait-il de même pour la population ? S'agirait-il exactement de la même relation linéaire ? Pourrions-nous prédire la réponse d'un nouvel individu qui ne fait pas partie de notre échantillon ?

Ces questions d'inférence et de prédiction se posent si nous pensons qu'un diagramme de dispersion reflète la relation sous-jacente entre les deux variables représentées, mais qu'il ne la spécifie pas complètement. Par exemple, un nuage de points représentant le poids à la naissance en fonction du nombre de jours de gestation nous montre la relation précise entre les deux variables dans notre échantillon ; mais nous pouvons nous demander si cette relation est vraie, ou presque, pour tous les bébés de la population dont l'échantillon a été tiré, ou même pour les bébés en général.

Comme toujours, la réflexion déductive commence par un examen minutieux des hypothèses relatives aux données. Les ensembles d'hypothèses sont appelés modèles. Les ensembles d'hypothèses sur le caractère aléatoire des diagrammes de dispersion à peu près linéaires sont appelés modèles de régression.

Un modèle de régression

En résumé, ces modèles affirment que la relation sous-jacente entre les deux variables est parfaitement linéaire ; cette ligne droite est le signal que nous aimerions identifier. Cependant, nous ne sommes pas en mesure de voir clairement la ligne. Ce que nous voyons, ce sont des points dispersés autour de la ligne. En chacun de ces points, le signal a été contaminé par du bruit aléatoire. Notre objectif déductif est donc de séparer le signal du bruit.

Plus précisément, le modèle de régression spécifie que les points du nuage de points sont générés au hasard comme suit.

• La relation entre $x$ et $y$ est parfaitement linéaire. Nous ne pouvons pas voir cette "vraie ligne", mais elle existe.
• Le diagramme de dispersion est créé en prenant des points sur la ligne et en les poussant verticalement hors de la ligne, soit au-dessus, soit au-dessous, comme suit :
  • Pour chaque $x$, trouver le point correspondant sur la vraie ligne (c'est le signal), puis générer le bruit ou l'erreur.
  • Les erreurs sont tirées au hasard avec remplacement à partir d'une population d'erreurs qui a une distribution normale avec une moyenne de 0.
  • Créez un point dont la coordonnée horizontale est $x$ et dont la coordonnée verticale est "la hauteur de la vraie ligne à $x$, plus l'erreur".
• Enfin, effacez la vraie droite de la dispersion et affichez uniquement les points créés.

def standard_units(any_numbers):
    "Convert any array of numbers to standard units."
    return (any_numbers - np.mean(any_numbers)) / np.std(any_numbers)

def correlation(t, x, y):
    return np.mean(standard_units(t[x]) * standard_units(t[y]))

def slope(table, x, y):
    r = correlation(table, x, y)
    return r * np.std(table[y]) / np.std(table[x])

def intercept(table, x, y):
    a = slope(table, x, y)
    return np.mean(table[y]) - a * np.mean(table[x])

def fit(table, x, y):
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

def scatter_fit(table, x, y):
    plots.scatter(table[x], table[y], s=20)
    plots.plot(table[x], fit(table, x, y), lw=2, color='gold')
    plots.xlabel(x)
    plots.ylabel(y)

Sur la base de ce diagramme de dispersion, comment devrions-nous estimer la vraie ligne ? La meilleure ligne que l'on puisse tracer dans un diagramme de dispersion est la ligne de régression. La droite de régression est donc une estimation naturelle de la vraie droite.

La simulation ci-dessous montre à quel point la droite de régression est proche de la vraie droite. Le premier panneau montre comment le diagramme de dispersion est généré à partir de la vraie droite. Le deuxième montre le diagramme de dispersion que nous voyons. Le troisième montre la ligne de régression à travers le diagramme. Le quatrième montre à la fois la droite de régression et la droite vraie.

Pour exécuter la simulation, appelez la fonction draw_and_compare avec trois arguments : la pente de la vraie droite, l'ordonnée à l'origine de la vraie droite et la taille de l'échantillon.

Exécutez la simulation plusieurs fois, avec différentes valeurs pour la pente et l'ordonnée à l'origine de la vraie droite, et différentes tailles d'échantillon. Étant donné que tous les points sont générés conformément au modèle, vous constaterez que la droite de régression est une bonne estimation de la vraie droite si la taille de l'échantillon est modérément grande.

def draw_and_compare(true_slope, true_int, sample_size):
    x = np.random.normal(50, 5, sample_size)
    xlims = np.array([np.min(x), np.max(x)])
    eps = np.random.normal(0, 6, sample_size)
    y = (true_slope * x + true_int) + eps
    tyche = pd.DataFrame({
        'x': x,
        'y': y
    })

    plots.figure(figsize=(6, 16))
    
    # Plot 1: True line and generated points
    plots.subplot(4, 1, 1)
    plots.scatter(tyche['x'], tyche['y'], s=20)
    plots.plot(xlims, true_slope * xlims + true_int, lw=2, color='green')
    plots.title('True Line, and Points Created')

    # Plot 2: Points only
    plots.subplot(4, 1, 2)
    plots.scatter(tyche['x'], tyche['y'], s=20)
    plots.title('What We Get to See')

    # Plot 3: Regression line
    plots.subplot(4, 1, 3)
    scatter_fit(tyche, 'x', 'y')
    plots.xlabel("")
    plots.ylabel("")
    plots.title('Regression Line: Estimate of True Line')

    # Plot 4: Regression line with true line
    plots.subplot(4, 1, 4)
    scatter_fit(tyche, 'x', 'y')
    plots.ylabel("")
    plots.plot(xlims, true_slope * xlims + true_int, lw=2, color='green')
    plots.title("Regression Line and True Line")

# The true line,
# the points created,
# and our estimate of the true line.
# Arguments: true slope, true intercept, number of points

draw_and_compare(4, -5, 10)

Dans la réalité, bien sûr, nous ne verrons jamais la vraie ligne. La simulation montre que si le modèle de régression semble plausible et si nous disposons d'un grand échantillon, la ligne de régression est une bonne approximation de la ligne réelle.

Inférence pour la pente réelle

Nos simulations montrent que si le modèle de régression est valable et que la taille de l'échantillon est importante, la droite de régression sera probablement proche de la vraie droite. Cela nous permet d'estimer la pente de la vraie droite.

Nous utiliserons notre échantillon familier de mères et de leurs nouveau-nés pour développer une méthode d'estimation de la pente de la vraie droite. Tout d'abord, voyons si nous pensons que le modèle de régression est un ensemble d'hypothèses approprié pour décrire la relation entre le poids de naissance et le nombre de jours de gestation.

def standard_units(any_numbers):
    "Convert any array of numbers to standard units."
    return (any_numbers - np.mean(any_numbers)) / np.std(any_numbers)

def correlation(t, x, y):
    return np.mean(standard_units(t[x]) * standard_units(t[y]))

def slope(table, x, y):
    r = correlation(table, x, y)
    return r * np.std(table[y]) / np.std(table[x])

def intercept(table, x, y):
    a = slope(table, x, y)
    return np.mean(table[y]) - a * np.mean(table[x])

def fit(table, x, y):
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

def residual(table, x, y):
    return table[y] - fit(table, x, y)

def scatter_fit(table, x, y):
    plots.scatter(table[x], table[y], s=20)
    plots.plot(table[x], fit(table, x, y), lw=2, color='gold')
    plots.xlabel(x)
    plots.ylabel(y)

def draw_and_compare(true_slope, true_int, sample_size):
    x = np.random.normal(50, 5, sample_size)
    xlims = np.array([np.min(x), np.max(x)])
    eps = np.random.normal(0, 6, sample_size)
    y = (true_slope * x + true_int) + eps
    tyche = pd.DataFrame({
        'x': x,
        'y': y
    })

    plots.figure(figsize=(6, 16))
    
    # Plot 1: True line and points created
    plots.subplot(4, 1, 1)
    plots.scatter(tyche['x'], tyche['y'], s=20)
    plots.plot(xlims, true_slope * xlims + true_int, lw=2, color='green')
    plots.title('True Line, and Points Created')

    # Plot 2: Points only
    plots.subplot(4, 1, 2)
    plots.scatter(tyche['x'], tyche['y'], s=20)
    plots.title('What We Get to See')

    # Plot 3: Regression line
    plots.subplot(4, 1, 3)
    scatter_fit(tyche, 'x', 'y')
    plots.xlabel("")
    plots.ylabel("")
    plots.title('Regression Line: Estimate of True Line')

    # Plot 4: Regression line and true line
    plots.subplot(4, 1, 4)
    scatter_fit(tyche, 'x', 'y')
    plots.ylabel("")
    xlims = np.array([np.min(tyche['x']), np.max(tyche['x'])])
    plots.plot(xlims, true_slope * xlims + true_int, lw=2, color='green')
    plots.title("Regression Line and True Line")

baby = pd.read_csv(path_data + 'baby.csv')

scatter_fit(baby, 'Gestational Days', 'Birth Weight')

correlation(baby, 'Gestational Days', 'Birth Weight')

0.40754279338885196

Dans l'ensemble, la dispersion semble assez uniforme autour de la ligne, bien que certains points soient dispersés à la périphérie du nuage principal. La corrélation est de 0,4 et la ligne de régression a une pente positive.

Cela reflète-t-il le fait que la vraie ligne a une pente positive ? Pour répondre à cette question, voyons si nous pouvons estimer la véritable pente. Nous disposons certainement d'une estimation de celle-ci : la pente de notre ligne de régression. Cela représente environ 0,47 once par jour.

slope(baby, 'Gestational Days', 'Birth Weight')

0.4665568769492164

Mais si le diagramme de dispersion avait été différent, la ligne de régression aurait été différente et aurait pu avoir une pente différente. Comment déterminer à quel point la pente aurait pu être différente ?

Nous avons besoin d'un autre échantillon de points, afin de pouvoir tracer la ligne de régression à travers le nouveau nuage de points et trouver sa pente. Mais où trouver un autre échantillon ?

Vous l'avez deviné, nous allons bootstrapper notre échantillon original. Nous obtiendrons ainsi un diagramme de dispersion bootstrap, à travers lequel nous pourrons tracer une ligne de régression.

Bootstrap du nuage de points

Nous pouvons simuler de nouveaux échantillons en procédant à un échantillonnage aléatoire avec remplacement à partir de l'échantillon original, autant de fois que la taille de l'échantillon original. Chacun de ces nouveaux échantillons nous donnera un diagramme de dispersion. Nous appellerons cela un diagramme de dispersion bootstrapé, et pour faire court, nous appellerons l'ensemble du processus bootstrapping the scatter plot (littéralement : "l'amorçage du diagramme de dispersion").

Voici le diagramme de dispersion original de l'échantillon et quatre répétitions de la procédure de rééchantillonnage bootstrap. Remarquez que les diagrammes de dispersion rééchantillonnés sont en général un peu plus clairsemés que l'original. Cela s'explique par le fait que certains des points originaux ne sont pas sélectionnés dans les échantillons.

plots.figure(figsize=(8, 18))

# Plot the original sample
plots.subplot(5, 1, 1)
plots.scatter(baby['Gestational Days'], baby['Birth Weight'], s=10)
plots.xlim([150, 400])
plots.title('Original sample')

# Plot 4 bootstrap samples
for i in np.arange(1, 5, 1):
    plots.subplot(5, 1, i + 1)
    rep = baby.sample(frac=1, replace=True)  # Bootstrap sample
    plots.scatter(rep['Gestational Days'], rep['Birth Weight'], s=10)
    plots.xlim([150, 400])
    plots.title('Bootstrap sample ' + str(i))

Estimation de la pente réelle

Nous pouvons extraire le diagramme de dispersion un grand nombre de fois et tracer une ligne de régression à travers chaque diagramme extrapolé. Chacune de ces lignes a une pente. Nous pouvons simplement collecter toutes les pentes et dessiner leur histogramme empirique. Rappelons que par défaut, la méthode sample tire au hasard avec remplacement, le même nombre de fois que le nombre de lignes dans le tableau. En d'autres termes, sample génère un échantillon bootstrap par défaut.

slopes = []

for i in np.arange(5000):
    bootstrap_sample = baby.sample(frac=1, replace=True)  # Bootstrap sample
    bootstrap_slope = slope(bootstrap_sample, 'Gestational Days', 'Birth Weight')
    slopes.append(bootstrap_slope)

# Convert slopes to a pandas DataFrame for visualization
slopes_df = pd.DataFrame({'Bootstrap Slopes': slopes})

# Plot histogram of bootstrap slopes
plots.hist(slopes_df['Bootstrap Slopes'], bins=20, edgecolor='black')
plots.xlabel('Bootstrap Slopes')
plots.ylabel('Frequency')
plots.title('Distribution of Bootstrap Slopes')

Text(0.5, 1.0, 'Distribution of Bootstrap Slopes')

Nous pouvons alors construire un intervalle de confiance approximatif de 95 % pour la pente de la vraie ligne, en utilisant la méthode du percentile bootstrap. L'intervalle de confiance s'étend du 2,5ème percentile au 97,5ème percentile des 5000 pentes bootstrapées.

left = np.percentile(slopes, 2.5)
right = np.percentile(slopes, 97.5)
left, right

(0.38409331573673783, 0.5570125075185108)

Un intervalle de confiance approximatif de 95 % pour la pente réelle s'étend d'environ 0,38 once par jour à environ 0,56 once par jour.

Une fonction pour amorcer la pente

Rassemblons toutes les étapes de notre méthode d'estimation de la pente et définissons une fonction bootstrap_slope qui les exécute. Ses arguments sont le nom du tableau et les étiquettes des variables de prédiction et de réponse, ainsi que le nombre souhaité de réplications bootstrap. Dans chaque réplication, la fonction bootstrape le nuage de points original et calcule la pente de la droite de régression résultante. Elle dessine ensuite l'histogramme de toutes les pentes générées et imprime l'intervalle composé des "95 % du milieu" des pentes.

def bootstrap_slope(table, x, y, repetitions):
    # Initialize an empty list to store bootstrap slopes
    slopes = []

    # For each repetition:
    for i in range(repetitions):
        # Bootstrap the sample
        bootstrap_sample = table.sample(frac=1, replace=True)
        # Calculate the slope of the regression line
        bootstrap_slope = slope(bootstrap_sample, x, y)
        # Append to the list of slopes
        slopes.append(bootstrap_slope)

    # Convert slopes to a NumPy array for percentile calculations
    slopes = np.array(slopes)

    # Find the endpoints of the 95% confidence interval for the true slope
    left = np.percentile(slopes, 2.5)
    right = np.percentile(slopes, 97.5)

    # Slope of the regression line from the original sample
    observed_slope = slope(table, x, y)

    # Plot the histogram of bootstrap slopes
    plots.hist(slopes, bins=20, edgecolor='black')
    plots.axvline(x=left, color='yellow', lw=2, label='2.5% CI')
    plots.axvline(x=right, color='yellow', lw=2, label='97.5% CI')
    plots.xlabel('Bootstrap Slopes')
    plots.ylabel('Frequency')
    plots.title('Distribution of Bootstrap Slopes')
    plots.legend()

    # Display results
    print('Slope of regression line:', observed_slope)
    print('Approximate 95%-confidence interval for the true slope:')
    print(left, right)

Lorsque nous appelons bootstrap_slope pour trouver un intervalle de confiance pour la vraie pente lorsque la variable réponse est le poids de naissance et le prédicteur est le nombre de jours de gestation, nous obtenons un intervalle très proche de celui que nous avons obtenu plus tôt : environ 0,38 onces par jour à 0,56 onces par jour.

bootstrap_slope(baby, 'Gestational Days', 'Birth Weight', 5000)

Slope of regression line: 0.4665568769492164
Approximate 95%-confidence interval for the true slope:
0.3804724830243045 0.5575795542762265

Maintenant que nous disposons d'une fonction qui automatise notre processus d'estimation de la pente de la vraie ligne dans un modèle de régression, nous pouvons également l'utiliser pour d'autres variables.

Par exemple, examinons la relation entre le poids à la naissance et la taille de la mère. Les femmes plus grandes ont-elles tendance à avoir des bébés plus lourds ?

Le modèle de régression semble raisonnable, d'après le diagramme de dispersion, mais la corrélation n'est pas élevée. Elle n'est que d'environ 0,2.

scatter_fit(baby, 'Maternal Height', 'Birth Weight')

correlation(baby, 'Maternal Height', 'Birth Weight')

0.20370417718968062

Comme précédemment, nous pouvons utiliser bootstrap_slope pour estimer la pente de la vraie ligne dans le modèle de régression.

bootstrap_slope(baby, 'Maternal Height', 'Birth Weight', 5000)

Slope of regression line: 1.4780193519284357
Approximate 95%-confidence interval for the true slope:
1.0406077026402254 1.9088665289902151

Un intervalle de confiance de 95 % pour la pente réelle s'étend d'environ 1 once par pouce à environ 1,9 once par pouce.

La pente réelle pourrait-elle être de 0 ?

Supposons que nous pensions que nos données suivent le modèle de régression et que nous ajustions la droite de régression pour estimer la vraie droite. Si la ligne de régression n'est pas parfaitement plate, comme c'est presque toujours le cas, nous observerons une certaine association linéaire dans le nuage de points.

Mais que se passe-t-il si cette observation est fallacieuse ? En d'autres termes, que se passerait-il si la véritable ligne était plate - c'est-à-dire s'il n'y avait pas de relation linéaire entre les deux variables - et si l'association que nous avons observée était simplement due au hasard dans la génération des points qui forment notre échantillon ?

Voici une simulation qui illustre la raison pour laquelle cette question se pose. Nous allons à nouveau appeler la fonction draw_and_compare, en exigeant cette fois que la vraie ligne ait une pente de 0. Notre objectif est de voir si notre droite de régression présente une pente différente de 0.

Rappelez-vous que les arguments de la fonction draw_and_compare sont la pente et l'ordonnée à l'origine de la vraie droite, ainsi que le nombre de points à générer.

draw_and_compare(0, 10, 25)

Effectuez la simulation plusieurs fois, en maintenant la pente de la ligne réelle à 0 à chaque fois. Vous remarquerez que si la pente de la vraie droite est égale à 0, la pente de la droite de régression n'est généralement pas égale à 0. La ligne de régression est tantôt ascendante, tantôt descendante, ce qui nous donne à chaque fois l'impression erronée que les deux variables sont corrélées.

Afin de déterminer si la pente que nous observons est réelle ou non, nous souhaitons tester les hypothèses suivantes :

**Hypothèse nulle : la pente de la vraie droite est 0.

**Hypothèse alternative : la pente de la vraie droite est différente de 0.

Nous sommes bien placés pour le faire. Puisque nous pouvons construire un intervalle de confiance à 95 % pour la pente réelle, il nous suffit de vérifier si l'intervalle contient 0.

Si ce n'est pas le cas, nous pouvons rejeter l'hypothèse nulle (avec le seuil de 5 % pour la valeur P).

Si l'intervalle de confiance pour la vraie pente contient 0, nous n'avons pas assez de preuves pour rejeter l'hypothèse nulle. Il se peut que la pente que nous observons soit fausse.

Utilisons cette méthode dans un exemple. Supposons que nous essayons d'estimer le poids du bébé à la naissance en fonction de l'âge de la mère. Sur la base de l'échantillon, la pente de la ligne de régression pour l'estimation du poids de naissance en fonction de l'âge de la mère est positive, d'environ 0,08 once par an.

slope(baby, 'Maternal Age', 'Birth Weight')

0.08500766941582515

Bien que la pente soit positive, elle est assez faible. La ligne de régression est si proche de l'horizontale que l'on peut se demander si la vraie ligne est horizontale.

scatter_fit(baby, 'Maternal Age', 'Birth Weight')

Nous pouvons utiliser bootstrap_slope pour estimer la pente de la vraie ligne. Le calcul montre qu'un intervalle de confiance bootstrap approximatif de 95 % pour la vraie pente a une extrémité gauche négative et une extrémité droite positive - en d'autres termes, l'intervalle contient 0.

bootstrap_slope(baby, 'Maternal Age', 'Birth Weight', 5000)

Slope of regression line: 0.08500766941582515
Approximate 95%-confidence interval for the true slope:
-0.10200994728774188 0.271715869022962

Comme l'intervalle contient 0, nous ne pouvons pas rejeter l'hypothèse nulle selon laquelle la pente de la véritable relation linéaire entre l'âge maternel et le poids du bébé à la naissance est 0. Sur la base de cette analyse, il ne serait pas judicieux de prédire le poids à la naissance sur la base du modèle de régression avec l'âge maternel comme variable prédictive.

Intervalles de prédiction

L'une des principales utilisations de la régression consiste à faire des prédictions pour un nouvel individu qui ne faisait pas partie de notre échantillon initial mais qui est similaire aux individus échantillonnés. Dans le langage du modèle, nous voulons estimer $y$ pour une nouvelle valeur de $x$.

Notre estimation est la hauteur de la vraie ligne à $x$. Bien entendu, nous ne connaissons pas la vraie ligne. Ce que nous avons comme substitut est la ligne de régression à travers notre échantillon de points.

La valeur ajustée à une valeur donnée de $x$ est l'estimation de régression de $y$ basée sur cette valeur de $x$. En d'autres termes, la valeur ajustée à une valeur donnée de $x$ est la hauteur de la droite de régression à cette valeur de $x$.

baby = pd.read_csv(path_data + 'baby.csv')

def standard_units(any_numbers):
    "Convert any array of numbers to standard units."
    return (any_numbers - np.mean(any_numbers)) / np.std(any_numbers)

def correlation(t, x, y):
    return np.mean(standard_units(t[x]) * standard_units(t[y]))

def slope(table, x, y):
    r = correlation(table, x, y)
    return r * np.std(table[y]) / np.std(table[x])

def intercept(table, x, y):
    a = slope(table, x, y)
    return np.mean(table[y]) - a * np.mean(table[x])

def fit(table, x, y):
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * table[x] + b

def residual(table, x, y):
    return table[y] - fit(table, x, y)

def scatter_fit(table, x, y):
    plots.scatter(table[x], table[y], s=20)
    plots.plot(table[x], fit(table, x, y), lw=2, color='gold')
    plots.xlabel(x)
    plots.ylabel(y)

Supposons que nous essayons de prédire le poids d'un bébé à la naissance en fonction du nombre de jours de gestation. Comme nous l'avons vu dans la section précédente, les données correspondent assez bien au modèle de régression et un intervalle de confiance de 95 % pour la pente de la vraie ligne ne contient pas 0. Il semble donc raisonnable de réaliser notre prédiction.

La figure ci-dessous montre où se situe la prédiction sur la ligne de régression. La ligne rouge se situe à $x = 300$.

scatter_fit(baby, 'Gestational Days', 'Birth Weight')

s = slope(baby, 'Gestational Days', 'Birth Weight')
i = intercept(baby, 'Gestational Days', 'Birth Weight')
fit_300 = s * 300 + i

# Scatter the point (300, fit_300) in red
plots.scatter(300, fit_300, color='red', s=20)

# Plot a vertical red line from 0 to fit_300
plots.plot([300, 300], [0, fit_300], color='red', lw=2)

# Set the y-axis limits
plots.ylim([0, 200])

(0.0, 200.0)

La hauteur du point où la ligne rouge touche la ligne de régression est la valeur ajustée à 300 jours de gestation.

La fonction fitted_value calcule cette hauteur. Comme les fonctions correlation, slope et intercept, ses arguments incluent le nom du tableau et les étiquettes des colonnes $x$ et $y$. Mais elle requiert également un quatrième argument, qui est la valeur de $x$ à laquelle l'estimation sera faite.

def fitted_value(table, x, y, given_x):
    a = slope(table, x, y)
    b = intercept(table, x, y)
    return a * given_x + b

La valeur ajustée à 300 jours de gestation est d'environ 129,2 onces. En d'autres termes, pour une grossesse d'une durée de 300 jours gestationnels, notre estimation du poids du bébé est d'environ 129,2 onces.

fit_300 = fitted_value(baby, 'Gestational Days', 'Birth Weight', 300)
fit_300

129.21292417031435

La variabilité de la prédiction

Nous avons mis au point une méthode permettant de prédire le poids de naissance d'un nouveau-né en fonction du nombre de jours de gestation, à l'aide des données de notre échantillon. Mais en tant que scientifiques des données, nous savons que l'échantillon aurait pu être différent. Si l'échantillon avait été différent, la ligne de régression aurait également été différente, tout comme notre prédiction. Pour déterminer la qualité de notre prédiction, nous devons nous faire une idée de la variabilité de la prédiction.

Pour ce faire, nous devons générer de nouveaux échantillons. Nous pouvons le faire en bootstrapant le diagramme de dispersion comme dans la section précédente. Nous ajusterons ensuite la ligne de régression au diagramme de dispersion dans chaque réplication et ferons une prédiction basée sur chaque ligne. La figure ci-dessous montre 10 lignes de ce type et le poids de naissance prédit correspondant à 300 jours de gestation.

x = 300

lines = pd.DataFrame(columns=['slope', 'intercept'])

for i in range(10):
    # Bootstrap sample
    rep = baby.sample(frac=1, replace=True)
    a = slope(rep, 'Gestational Days', 'Birth Weight')
    b = intercept(rep, 'Gestational Days', 'Birth Weight')
    lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)

# Calculate predictions at x
lines['prediction at x=' + str(x)] = lines['slope'] * x + lines['intercept']

# Define limits for x and corresponding y values
xlims = np.array([291, 309])
left = xlims[0] * lines['slope'] + lines['intercept']
right = xlims[1] * lines['slope'] + lines['intercept']
fit_x = x * lines['slope'] + lines['intercept']

# Plot each line and the predictions at x
for i in range(10):
    plots.plot(xlims, np.array([left[i], right[i]]), lw=1)
    plots.scatter(x, fit_x[i], s=30)

C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\1573641999.py:10: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)

Les prédictions varient d'une ligne à l'autre. Le tableau ci-dessous indique la pente et l'ordonnée à l'origine de chacune des 10 droites, ainsi que la prédiction.

lines

	slope	intercept	prediction at x=300
0	0.457480	-8.837608	128.406412
1	0.398475	8.215161	127.757613
2	0.508144	-22.835515	129.607803
3	0.431348	-0.881831	128.522600
4	0.452061	-7.190722	128.427449
5	0.394275	9.442054	127.724465
6	0.527881	-26.809522	131.554797
7	0.420601	0.599649	126.779902
8	0.504815	-21.705768	129.738683
9	0.536698	-30.433994	130.575268

Intervalle de prédiction Bootstrap

Si nous augmentons le nombre de répétitions du processus de rééchantillonnage, nous pouvons générer un histogramme empirique des prédictions. Cela nous permettra de créer un intervalle de prédictions, en utilisant la même méthode de percentile que celle utilisée pour créer un intervalle de confiance bootstrap pour la pente.

Définissons une fonction appelée bootstrap_prediction pour faire cela. La fonction prend cinq arguments :

• le nom du tableau
• les étiquettes des colonnes des variables prédicteur et réponse, dans cet ordre
• la valeur de $x$ à partir de laquelle la prédiction doit être effectuée
• le nombre désiré de répétitions bootstrap

Lors de chaque répétition, la fonction bootstrap reprend le nuage de points original et trouve la valeur prédite de $y$ sur la base de la valeur spécifiée de $x$. Plus précisément, elle appelle la fonction fitted_value que nous avons définie plus tôt dans cette section pour trouver la valeur ajustée à la valeur $x$ spécifiée.

Enfin, il dessine l'histogramme empirique de toutes les valeurs prédites, et imprime l'intervalle composé des "95% du milieu" des valeurs prédites. Il imprime également la valeur prédite basée sur la ligne de régression à travers le diagramme de dispersion original.

def bootstrap_prediction(table, x, y, new_x, repetitions):
    # Initialize an empty list to store predictions
    predictions = []

    # Bootstrap procedure
    for i in np.arange(repetitions):
        # Create a bootstrap sample
        bootstrap_sample = table.sample(frac=1, replace=True)
        # Get the regression prediction at new_x
        bootstrap_prediction = fitted_value(bootstrap_sample, x, y, new_x)
        # Append the prediction to the list
        predictions.append(bootstrap_prediction)
        
    # Convert predictions to a NumPy array for calculations
    predictions = np.array(predictions)

    # Find the ends of the approximate 95% prediction interval
    left = np.percentile(predictions, 2.5)
    right = np.percentile(predictions, 97.5)

    # Prediction based on the original sample
    original = fitted_value(table, x, y, new_x)

    # Display results: Plot the histogram of predictions
    plots.hist(predictions, bins=20, edgecolor='black')
    plots.xlabel('Predictions at x=' + str(new_x))
    plots.axvline(x=left, color='yellow', lw=4, label='2.5% CI')
    plots.axvline(x=right, color='yellow', lw=4, label='97.5% CI')
    plots.legend()

    # Print the results
    print('Height of regression line at x=' + str(new_x) + ':', original)
    print('Approximate 95%-confidence interval:')
    print(left, right)

bootstrap_prediction(baby, 'Gestational Days', 'Birth Weight', 300, 5000)

Height of regression line at x=300: 129.21292417031435
Approximate 95%-confidence interval:
127.28604701356356 131.2830905699535

La figure ci-dessus montre un histogramme empirique bootstrap du poids de naissance prédit d'un bébé à 300 jours de gestation, basé sur 5 000 répétitions du processus bootstrap. La distribution empirique est à peu près normale.

Un intervalle de prédiction approximatif de 95 % des scores a été construit en prenant le "milieu de 95 %" des prédictions, c'est-à-dire l'intervalle allant du 2,5ème centile au 97,5ème centile des prédictions. L'intervalle va d'environ 127 à environ 131. La prédiction basée sur l'échantillon original était d'environ 129, ce qui est proche du centre de l'intervalle.

L'effet de la modification de la valeur du prédicteur

La figure ci-dessous montre l'histogramme de 5 000 prédictions bootstrap à 285 jours de gestation. La prédiction basée sur l'échantillon original est d'environ 122 onces, et l'intervalle s'étend d'environ 121 onces à environ 123 onces.

bootstrap_prediction(baby, 'Gestational Days', 'Birth Weight', 285, 5000)

Height of regression line at x=285: 122.2145710160761
Approximate 95%-confidence interval:
121.15639820273395 123.28390764523947

Remarquez que cet intervalle est plus étroit que l'intervalle de prédiction à 300 jours de gestation. Cherchons-en la raison.

Le nombre moyen de jours de gestation est d'environ 279 jours :

np.mean(baby['Gestational Days'])

279.1013628620102

Ainsi, 285 est plus proche du centre de la distribution que 300. En règle générale, les lignes de régression basées sur les échantillons bootstrap sont plus proches les unes des autres près du centre de la distribution de la variable prédictive. Par conséquent, toutes les valeurs prédites sont également plus proches les unes des autres. Cela explique la largeur réduite de l'intervalle de prédiction.

Vous pouvez le constater dans la figure ci-dessous, qui montre les prédictions à $x = 285$ et $x = 300$ pour chacune des dix réplications bootstrap. En règle générale, les lignes sont plus éloignées à $x = 300$ qu'à $x = 285$, et les prédictions à $x = 300$ sont donc plus variables.

x1 = 300
x2 = 285

lines = pd.DataFrame(columns=['slope', 'intercept'])

# Bootstrap process to generate slopes and intercepts
for i in range(10):
    rep = baby.sample(frac=1, replace=True)
    a = slope(rep, 'Gestational Days', 'Birth Weight')
    b = intercept(rep, 'Gestational Days', 'Birth Weight')
    lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)

# Define x-limits
xlims = np.array([260, 310])

# Calculate the left and right points for each line
left = xlims[0] * lines['slope'] + lines['intercept']
right = xlims[1] * lines['slope'] + lines['intercept']

# Calculate predictions for x1 and x2
fit_x1 = x1 * lines['slope'] + lines['intercept']
fit_x2 = x2 * lines['slope'] + lines['intercept']

# Plot all lines and predictions
plots.xlim(xlims)
for i in range(10):
    plots.plot(xlims, np.array([left[i], right[i]]), lw=1)
    plots.scatter(x1, fit_x1[i], s=30)
    plots.scatter(x2, fit_x2[i], s=30)

C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)
C:\Users\guenser\AppData\Local\Temp\ipykernel_12464\746407846.py:11: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  lines = lines.append({'slope': a, 'intercept': b}, ignore_index=True)

Mise en garde

Toutes les prédictions et tous les tests que nous avons effectués dans ce chapitre supposent que le modèle de régression est valable. Plus précisément, les méthodes supposent que le nuage de points ressemble à des points générés en commençant par des points situés sur une ligne droite, puis en les poussant hors de la ligne par l'ajout d'un bruit normal aléatoire.

Si le nuage de points ne ressemble pas à cela, il se peut que le modèle ne corresponde pas aux données. Si le modèle ne tient pas, les calculs qui supposent que le modèle est vrai ne sont pas valables.

Par conséquent, nous devons d'abord décider si le modèle de régression s'applique à nos données, avant de commencer à faire des prédictions basées sur le modèle ou de tester des hypothèses sur les paramètres du modèle. Une méthode simple consiste à faire ce que nous avons fait dans cette section, c'est-à-dire dessiner le diagramme de dispersion des deux variables et voir s'il est à peu près linéaire et uniformément réparti autour d'une ligne. Nous devrions également exécuter les diagnostics que nous avons développés dans la section précédente en utilisant le diagramme des résidus.

Crédits

Ce cours est inspiré du cours data8 donné à UC Berkeley et en ré-utilise avec certaines modifications une partie des matériels (ces matériels sont généreusement mis à disposition publiquement sous licence Creative Commons avec attribution, consultez https://www.data8.org pour plus d'informations.

La dernière partie s'inspire également du cours « Exploring and Modeling High-Dimensional Data » de « Carpentries Incubator » et en ré-utilise avec certaines modifications une partie des matériels (ces matériels sont généreusement mis à disposition publiquement sous licence Creative Commons avec attribution, consultez https://carpentries-incubator.github.io/high-dimensional-analysis-in-python pour plus d'informations.